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Abstract 


Artificially  generated  plasmas  may  be  employed  to  alter  the  propagation  characteristics 
of  electromagnetic  waves.  The  purpose  of  this  report  is  to  study  the  propagation  of 
electromagnetic  waves  in  an  electron  beam  generated  plasma.  To  understand  the  physics 
related  to  this  concept  requires  the  development  of  computational  tools  dealing  with  a 
plasma  created  by  an  electron  beam,  an  assessment  of  the  temporal  and  spatial  evolution 
of  the  plasma,  and  a  characterization  of  the  refraction  and  attenuation  of  electromagnetic 
(EM)  waves  in  a  collisional  plasma.  Three  computer  programs  were  developed  to 
characterize  the  effectiveness  of  an  electron  beam  generated  plasma  in  refracting  and 
attenuating  an  EM  wave.  The  spatial  extent  and  density  distribution  of  a  plasma 
generated  by  a  relativistic  electron  beam  were  determined  using  an  axisymmetric  Monte 
Carlo  model.  This  plasma  density  distribution  was  used  as  a  source  term  in  the  second 
code,  a  temporal  solution  of  the  plasma  evolution  based  on  a  time  dependent  analysis  of 
the  plasma  rate  equations.  The  third  code  developed,  evaluates  the  attenuation  and 
refraction  of  an  EM  wave  in  the  resulting  plasma  by  using  a  ray  tracing  method  based  on 
the  eikonal  approach  of  Sommerfeld.  The  theoretical  foundation  and  validation 
procedures  are  presented  for  each  program.  A  limited  exploration  of  the  dependence  of 
the  plasma  distribution  on  neutral  densities  and  the  electron  beam  energies  was 
performed.  For  neutral  densities  corresponding  to  5  km  altitude,  the  plasma  longitudinal 
extent  ranged  from  52  to  868  cm  and  the  radial  extent  ranged  from  18  to  292  cm  for 
initial  electron  energies  between  100  keV  and  1  MeV  respectively.  Plasma  chemistry 
plays  a  critical  role  in  determining  the  electron  plasma  density  and  dictates  the  beam 
format  required  to  achieve  a  desired  level  of  EM  wave  attenuation. 
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ATTENUATION  AND  REFRACTION  OF  AN  ELECTROMAGNETIC  WAVE 
IN  AN  ELECTRON  BEAM  GENERATED  PLASMA 

I.  Introduction 

The  purpose  of  this  study  is  to  examine  the  attenuation  and  refraction  of  an  EM 
wave  traversing  an  electron  beam  generated  plasma.  This  chapter  starts  with  a  more 
detailed  statement  of  the  objectives  of  the  study  and  then  describes  the  general  approach 
taken  to  achieve  those  objectives.  The  last  part  of  the  chapter  gives  background 
information  on  plasma  characteristics  that  will  be  useful  for  the  remainder  of  the  study. 

Chapter  II  first  summarizes  the  theory  of  the  refraction  and  attenuation  of  an  EM 
wave  propagating  through  a  collisional  plasma.  The  computer  program,  written  to 
evaluate  the  refraction  and  attenuation  of  an  EM  wave  traversing  an  electron  beam 
generated  plasma,  is  then  discussed. 

Chapter  III  provides  background  theory  on  electron  collision  cross  sections  that 
are  used  to  model  the  plasma  generation.  The  theoretical  cross  sections  are  compared 
with  experimental  values  compiled  by  the  National  Institute  of  Standards  and  Technology 
(NIST).  The  comparison  is  done  to  validate  the  theoretical  models  as  well  as  insure  that 
they  are  implemented  correctly. 

In  Chapter  IV,  the  simple  electron  beam  propagation  model  used  to  obtain 
bounding  values  for  the  densities  and  spatial  extent  of  the  plasma  is  introduced.  The 
Monte  Carlo  based  program  used  to  determine  the  plasma  density  and  spatial  distribution 
is  then  discussed.  The  plasma  loss  mechanisms  and  a  model  for  estimating  the  loss  in  the 
plasma  density  are  then  considered  and  addressed. 
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Chapter  V  presents  a  demonstration  of  the  capabilities  of  the  computer 
simulations  that  were  developed  in  the  first  four  chapters.  Results  from  the  computer 
simulations  are  discussed  and  conclusions  are  drawn.  Finally,  the  limitations  to  the  study 
and  recommendations  for  future  work  are  discussed. 

Objectives 

The  objective  of  this  study  is  to  examine  the  attenuation  and  refraction  of  an  EM 
wave  traversing  an  electron  beam  generated  plasma.  The  significance  of  the  spatially 
dependent  attenuation  will  be  cast  in  terms  of  a  spatially  averaged  attenuation  of  the 
incident  EM  wave. 

Specifically  this  study  was  designed  to  develop  tools  to  determine: 

1 .  The  electron  density  distribution  of  a  plasma  generated  by  a  relativistic  electron 
beam. 

2.  The  temporal  and  spatial  evolution  of  the  plasma  density  accounting  for  attachment 
and  recombination. 

3.  The  spatial  attenuation  and  refraction  of  an  EM  wave  with  finite  spatial  extent  due  to 
the  temporally  evolving  electron  density  distribution  resulting  from  a  relativistic 
electron  beam.  The  spatial  attenuation  and  refraction  is  determined  as  a  function  of 
certain  electron  beam  and  environmental  parameters  such  as  power,  initial  electron 
energy,  air  density  and  temperature. 

The  determination  of  the  electron  density  and  spatial  distribution  caused  by  a 
relativistic  electron  beam  ionizing  the  air  was  approached  in  two  phases.  The  first  phase 
bounded  the  problem  by  examining  the  forward  scattering  case.  This  case  assumes  that 
the  electrons  are  not  scattered  laterally  and  all  energy  lost  by  the  initial  electrons  results 
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in  the  formation  of  new  electrons  through  ionization.  These  assumptions  result  in  an 
electron  beam  that  travels  straight  through  the  air  and  decays  only  in  energy  as  it 
propagates.  The  decay  in  energy  is  due  to  the  ionizing  collisions  that  occur  when  the 
incident  electron,  an  electron  from  the  electron  beam,  impacts  the  electrons  of  a  neutral 
molecule.  For  these  simplified  cases,  ionizing  collisions  result  in  a  loss  of  energy  to  the 
incident  electron  and  the  formation  of  a  new  electron  that  has  either  zero  energy  or  half 
the  energy  of  the  incident  electron.  The  second  phase  accurately  determines  the  electron 
density  distribution  using  a  Monte  Carlo  simulation.  In  the  Monte  Carlo  simulation,  a 
triple  differential  cross-section  (TDCS)  developed  by  Mott  was  used  to  determine  the 
angular  scattering  of  the  incident  electron,  the  amount  of  energy  imparted  to  the  ejected 
electron,  and  the  ejection  angle  of  the  liberated  electron  after  the  collision.  The  results  of 
the  Monte  Carlo  simulation  are  smoothed  using  group  statistics  and  used  to  determine  a 
two-dimensional  electron  density  distribution  for  the  plasma. 

The  variation  in  plasma  density  over  time  was  modeled  from  a  rate  equation 
standpoint  using  differential  equations  developed  from  the  various  attachment, 
recombination,  and  detachment  processes  that  occur  in  the  plasma.  A  Runge-Kutta 
numerical  method  was  employed  to  solve  seventeen  first  order  non-linear  differential 
equations.  Those  equations  describe  the  temporal  evolution  of  the  concentrations  of  the 
various  species  in  the  plasma.  To  simplify  the  problem,  only  electron  densities  and 
densities  of  atomic  and  molecular  nitrogen  and  oxygen  and  their  respective  positive  and 
negative  ions  were  considered  in  the  calculations.  The  results  of  these  calculations  were 
used  to  modify  the  plasma  density,  so  that  it  reflected  the  loss  of  electrons  due  to 
attachment  and  recombination  processes. 
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The  attenuation  and  refraction  of  the  EM  wave  is  calculated  by  using  a  ray  tracing 
method  based  on  the  Sommerfeld  eikonal  method.  This  method  propagates  the  wave 
through  the  plasma  iteratively  and  determines  the  amount  of  refraction  the  wave 
undergoes  based  on  the  index  of  refraction  of  the  plasma  and  its  gradient.  The  amount 
the  EM  wave  is  attenuated,  in  general,  depends  on  the  frequency  of  the  EM  wave  and  the 
electron,  positive  and  negative  ion,  and  neutral  density  of  the  plasma  and  their  respective 
temperatures. 


Background 


Tonks  and  Langmuir  used  the  word  “plasma”  in  1929  “to  designate  that  portion  of 
an  arc -type  discharge  in  which  the  densities  of  ions  and  electrons  are  high  but 
substantially  equal”(Sturrock,  1994:6).  However,  the  term  plasma  has  been  broadened  to 
describe  the  fourth  state  of  matter  in  which  a  large  number  of  the  atoms  or  molecules  of  a 
gas  have  been  ionized  or  have  an  electrical  charge.  Plasma  also  has  the  characteristic  of 
being  quasi-neutral  and  exhibiting  collective  effects.  A  parameter  that  is  commonly  used 
to  describe  the  collective  effects  of  a  plasma  is  the  plasma  frequency.  The  plasma 
frequency  describes  the  maximum  undamped  frequency  at  which  the  electrons  oscillate 
in  a  plasma.  The  plasma  frequency  for  an  electron  and  positive  ion  plasma  is  described 
by  the  equation 


0) 


p 


Am  .e 


m„ 


(1) 


where 

cop  =  plasma  frequency 
ne  =  electron  density 
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To  evaluate  the  effects  of  a  plasma  on  EM  wave  propagation,  a  dispersion 
relationship  is  developed  for  the  plasma.  The  dispersion  relationship  enables  us  to 
determine  the  wavelength,  the  phase  and  group  velocity,  and  the  index  of  refraction  of  an 
EM  wave  in  a  plasma.  If  the  index  of  refraction  is  complex  then  the  EM  wave  will 
attenuate  as  it  traverses  the  plasma.  The  simplest  dispersion  relationship  is  associated 
with  a  collisionless,  cold  plasma  with  no  impressed  magnetic  field: 

a)1  =a)p2  +c2k2  (2) 

where 


c o  =  angular  frequency  of  the  EM  wave 
k  =  the  wave  number  {In  I  A) 

A  cold  plasma  is  a  plasma  in  which  the  thermal  velocities  of  the  constituents  are 
negligible  (Sturrock,  1994:73).  From  equation  (2),  we  obtain  an  expression  for  the  group 
velocity  of  an  EM  wave: 


v 


g 


dco  0) 

—  =  c\  1 - , 

dk  y  (o 


and  the  index  of  refraction  for  that  EM  wave  is 


n  = 


(3) 


(4) 


where 


vg  =  group  velocity 
n  =  index  of  refraction 

However,  plasmas  in  a  dense  gas,  such  as  air  at  atmospheric  pressure,  have  a  large 
collision  frequency  between  its  constituents,  therefore,  we  must  use  a  dispersion 
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relationship  for  a  collisional  plasma.  The  dispersion  relationship  for  a  collisional,  cold 
plasma  is  derived  using  an  effective  frictional  force  term  combined  with  the  forces  on  the 
electrons  due  to  an  EM  wave 


Fcoii  =~mvv 

where 


(5) 


v  =  collision  frequency 
v  =  particle  velocity 

Using  Maxwell  and  Lorentz’s  equations,  we  obtain  the  following  dispersion  relationship 


cor 


t  • v 
1  + 1  — 

1 0 


•  +  c  k 


2/2 


(6) 


Using  a  hard  sphere  approximation,  the  collision  frequency  of  electrons  with  neutral 
particles  is  given  by 


where 


(7) 


(s) 

V  nm 

v  =  collision  frequency  of  the  plasma 

g 

a  =  hard  sphere  radius  of  the  molecules  (a=1.2xl0‘  cm) 
v  =  average  thermal  velocity  of  the  electrons 
Nm  =  molecular  number  density  of  air 
T  =  temperature  of  the  electrons 

(Ginzburg,  1984:  41).  Equation  (6)  results  in  an  index  of  refraction  described  by 
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n  =  nr  +  int 


(9) 


where 


co„  -2od  „ (O  i/ 

n  =n+_p - 1 — )/4. 

r  ,2/,. ,2  > 


(0  ((0  +V  ) 


1  +  - 


2  4 


2  \2 


£W  (fi;  +  ir  ) 


ft;4  -2con2co2  y 

»,-=(!+  V  ,  p— r 


(ft;2+t>2)  ^ft;2 (rw2  +  w2  - ft;p2 )2  +  y2ft;/ 


(10) 


(11) 


Due  to  the  index  of  refraction  being  complex  in  a  collisional  plasma,  the  EM  wave  will 
attenuate  as  it  propagates  through  the  plasma  (Clemmow,  1976:188-189). 
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II.  EM  Wave  Propagation  in  a  Plasma 


Refraction  of  an  EM  Wave  in  a  Collisional  Plasma 
Snell’s  law  is: 

sin  01  =  n2  sin  02  (1) 

where 

nl  =  index  of  refraction  for  the  initial  medium 
n2  =  index  of  refraction  for  the  final  medium 
6X  =  angle  of  EM  wave  in  the  initial  medium  with  respect  to  normal 
02  =  angle  of  EM  wave  in  the  final  medium  with  respect  to  normal 
This  simple  equation  describes  the  refraction  of  an  EM  wave  as  it  passes  from  one 
medium  to  another  with  a  different  index  of  refraction  and  provides  the  foundation  of 
geometric  optics.  An  approximation  of  the  refraction  experienced  by  an  EM  wave  using 
Snell’s  law  can  be  obtained  by  considering  a  medium  that  slowly  varies  in  index  of 
refraction,  such  that  it  can  be  divided  into  discrete  layers.  The  refraction  of  the  EM  wave 
is  calculated  at  the  boundary  of  each  layer  resulting  in  a  curved  trajectory  of  the  EM 
wave  as  it  traverses  the  medium.  A  sample  trajectory  of  the  EM  wave  is  given  in  Figure 
1. 
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where 


na>nl>n1>  n3 


At  the  peak  of  the  ray  trajectory,  the  ray  is  parallel  to  the  layers  of  changing  index  of 
refraction  and  perpendicular  to  Vn .  At  the  top  of  the  trajectory,  Snell’s  law  fails  us 
because  it  indicates  that  the  ray  would  travel  parallel  to  the  layers  without  being  refracted 
because  the  index  of  refraction  is  no  longer  changing.  This  failure  is  due  to  Snell’s  law 
considering  the  EM  wave  to  be  a  ray  with  no  spatial  extent.  If  we  consider  the  spatial 

c 

extent  and  phase  velocity  of  the  EM  wave,  described  by  — ,  then  at  the  top  of  the 

n 

trajectory  the  lower  portion  of  the  EM  wave  will  travel  slower  than  the  upper  portion  of 
the  EM  wave,  hence  causing  the  wave  to  refract  downward  (see  Figure  2). 


Vn 


Figure  2.  Refraction  an  EM  Wave  Propagating  Perpendicular  to  Vn 

The  eikonal  method  (described  in  the  next  section),  unlike  Snell’s  law,  considers  the 
curvature  of  the  phase  front  of  an  EM  wave  as  it  propagates  through  an  inhomogeneous 
medium.  Therefore,  it  will  be  used  to  determine  the  trajectory  of  an  EM  wave  as  it 
travels  through  the  plasma. 

Rav  Tracing  using  the  Eikonal  Method 

Eikonal  is  the  name  given  to  a  function  that  describes  the  constant  phase  front  of  a 
wave.  The  most  commonly  used  eikonals  are  planar,  cylindrical,  spherical,  and 
quadratic.  Geometric  optics  primarily  utilizes  planar  waves  while  Fourier  optics  utilizes 
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planar,  spherical,  and  quadratic  phase  fronts  in  a  homogenous  medium  with 
discontinuities  (i.e.  lenses,  aperture  stops,  prisms,  etc.).  The  eikonal,  however,  is  not 
restricted  to  those  simple  geometric  forms  for  the  phase  front  and  in  an  inhomogeneous 
medium  it  may  become  very  complex.  The  next  section  describes  three  different 
methods  that  utilize  the  eikonal  function  to  determine  the  trajectory  of  a  wave  through  an 
inhomogeneous  medium. 


Comparison  of  Eikonal  Methods 

Three  different  methods  of  determining  the  refraction  of  a  wave  in  an 
inhomogeneous  medium  were  inspected  for  possible  use  as  a  means  of  propagating  an 
EM  wave  through  the  plasma.  All  three  of  the  methods  inspected  utilized  the  eikonal 
approach,  however,  each  ray  tracing  method  was  developed  differently.  The  following 
section  compares  these  propagation  methodologies. 

The  first  ray  tracing  method  evaluated  was  named  the  Sommerfeld  Iterative 
Method  (SIM)  because  it  was  based  on  the  curvature  vector  equation  developed  by 
Sommerfeld.  The  curvature  vector  of  a  ray  in  an  inhomogeneous  medium  is  described  by 


K  =  —  (?xVn)x? 
n 


(2) 


where 


K  =  curvature  vector 
s  =  ray  propagation  unit  vector 

(Sommerfeld,  1964:339)(See  Appendix  A  for  details  on  the  derivation  of  equation  (2)). 
The  magnitude  of  the  curvature  vector  is  the  radius  of  curvature  of  the  path  of  the  EM 
wave  as  it  is  refracted  in  the  medium.  In  Appendix  A,  it  is  shown  that  the  curl  of  the 
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eikonal  is  zero;  hence  the  change  in  the  eikonal  is  path  independent.  Using  that  fact  and 
the  relationship 

dL  =  RdO  (3) 

where 


R  =  1/ 


K 


(4) 


dL  =  vgdt 


(5) 


We  can  derive  a  first  order  differential  equation  that  describes  the  rate  of  change  of  the 
angle  of  the  ray  propagation  direction  in  the  lab  coordinates  in  an  inhomogeneous 
medium.  The  rate  of  change  of  the  angle,  6 ,  in  an  inhomogeneous  medium  is  given  by 


—  =  v5(x,y)|f(x,y)|  (6) 


and  the  rate  of  change  of  the  x  and  y  position  of  the  ray  is  determined  by  the  x  and  y 
components  of  the  group  velocity  using  the  differential  equations 


—  =  v  (x,y)Cos(0)  (7) 

dt  s 

=  v  (x,  y)Sin(6)  (8) 

dt  * 


Equations  (6),  (7),  and  (8)  are  then  used  to  trace  the  ray  path  in  an  inhomogeneous 
medium.  These  equations  appear  to  be  benign  at  first,  but  no  analytic  solution  has  been 
obtained  from  them  except  for  the  trivial  case  of  a  homogeneous  medium. 
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Figure  3.  Diagram  of  the  Sommerfeld  Radius  of  Curvature 
Haslegrove  also  developed  a  set  of  differential  equations  for  ray  tracing  in  an 
inhomogeneous  medium.  His  equations  are  derived  from  the  differential  form  of  Snell’s 
Law,  but  are  very  similar  to  the  equations  obtained  in  the  SIM.  The  three  first  order 
differential  equations  are 


dx  c  .  ...  _.  //1N  dn. 

—  =  —  ( nCos{0 )  +  Sin(0)—) 
dt  n  dd 


(9) 


4  =  4  (.nSinW-CostfA 
dt  n  da 


(10) 


£-4(fW>-£staW))  (id 

dt  n  dy  dx 

(Haslegrove,  1954:355-358).  One  notable  difference  between  the  two  ray  tracing 

c 

methods  is  that  the  propagation  velocity  is  the  phase  velocity  of  the  wave,  — ,  in 

n 


Haslegrove’s  equations,  where  in  the  SIM  the  group  velocity  from  Section  I,  equation  (4) 
is  used  for  the  propagation  velocity  of  the  wave.  Since,  the  group  velocity  is  the  rate  at 
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which  the  energy  and  information  travel  through  the  medium,  it  is  of  more  interest  to  us 
than  the  phase  velocity.  Equations  (9)  -  (11)  are  intended  for  use  in  a  curvilinear  plane,  a 


dn  c 

curved  two-dimensional  plane,  which  is  the  reason  for  the  7—  term.  If  the  term  — 

06  n 


is 


replaced  by  vg  and  the  —terms  are  set  to  zero  then  the  Sommerfeld  Iterative  Method 


and  the  Haslegrove  Method  are  identical. 

Budden  also  developed  an  analytic  expression  for  the  path  of  a  ray  in  a  linear  and 
exponentially  varying  plasma.  The  analytic  solution  presented  later  in  this  section  was 
obtained  by  using  an  integral  equation  Budden  developed  to  trace  the  path  of  a  ray  in  a 
medium  that  varies  in  index  of  refraction  in  only  one-dimension  (Budden,  1961:178). 
The  form  of  the  integral  is 


(i2) 

For  a  medium  varying  in  one  dimension  the  following  relationships  can  be  used 


S  =  nSin{6 ) 

(13) 

2  2  c 2 

^  -  fl 

(14) 

Therefore 

dq  _  S 
dS  q 


(15) 


x  = 


(16) 


For  a  medium  that  varies  linearly  in  z  and  an  EM  wave  that  has  angle  of  incidence,  6l ,  to 
the  medium,  the  value  of  q  becomes 
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(17) 


2  2  r\  Ofc 

q  —  cos  &i  ~J2 


where 


/=  frequency  of  the  EM  wave 

cop(z)2 

a  =  — - 

2  7IZ 

Substituting  (17)  into  (16)  and  integrating,  we  obtain  the  expression  for  the  path  of  the 
EM  wave  through  the  linearly  varying  plasma 

/ 2  sin  20,  2/ 2  sin  0,  I  ~  ~o& 

x  =  - - - Llcos20,  (18) 

a  a  v  f 


It  should  be  noted  that  the  Budden’s  equation  predicts  that  the  path  will  be  exactly 
parabolic.  Next  we  will  consider  a  plasma  which  exponentially  varies  in  density  in  one 
dimension  such  that  q  has  the  form 


q2  -  cos2  6,  - 


2nf 


taz 


(19) 


where 


a)p2=C0p2(0)eoz 

(20) 

1,  co2 
a-—  In — f — 

(21) 

z  coD  (0) 

Substituting  (19)  into  (16)  and  integrating  we  obtain  the  expression  for  the  path  of  the 
EM  wave  in  an  exponentially  varying  plasma 


x  = 


2sin#; 
a  cos  6, 


In 


1 

tan — 

2  2 


1 

tan — 0, 

2  1 


(22) 
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where 


sin  cp2  - 


<»,(0)  j<* 

— - - e2 

co  cos  (9; 


sin  ^  = 


^(0) 

CO  COS  0J 


(23) 


(24) 


SIM  and  Haslegrove’s  differential  equations  usually  do  not  result  in  analytic 
solutions,  but  can  be  used  to  determine  the  trajectory  of  the  ray  using  a  standard 
numerical  technique  for  solving  differential  equations.  Haslegrove’s  equations  are 
limited  to  use  in  a  one-dimensional,  curvilinear  plane  whereas  the  SIM  can  be  used  for  a 
two  or  even  three-dimensional  varying  plasma.  Budden’s  equation  for  the  limited  cases 
of  a  one-dimensional,  linearly  and  exponential  varying  plasma  results  in  an  analytic 
solution.  However,  for  more  complex  medium  an  analytic  solution  is  rarely  obtained. 
There  is  also  a  two  dimensional  version  of  Budden’s  integral  equations  (Budden 
1961:176),  which  can  be  used  to  determine  ray  trajectories  using  numerical  integration 
techniques.  This  method,  however,  provides  no  capabilities  above  what  has  already  been 
presented  in  this  section. 

Comparison  of  Analytic  and  Numeric  Ray  Tracing  Results 

To  validate  the  Sommerfeld  Iterative  Method,  a  comparison  of  trajectory  results 
was  performed  for  the  three  methods  described  in  the  chapter.  The  first  comparison  case 
examined  the  trajectory  of  an  EM  wave  in  a  plasma  linearly  varying  in  density  in  one- 
dimension.  Figure  4.a  shows  the  trajectory  of  an  EM  wave  for  all  three  eikonal  methods 
in  a  plasma  that  increases  in  density  linearly  with  increasing  y  values.  From  Figure  4.b, 
we  see  that  the  index  of  refraction  of  the  plasma  decreases  with  increasing  plasma 
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density.  We  expect  from  Snell’s  law  that  the  angle  of  the  ray  relative  to  the  y  axis  would 
increase  because  n2  is  less  than  nx  which  results  in  02  increasing  to  compensate. 

Therefore,  the  ray  refracts  in  the  direction  of  Vn  over  the  entire  path  of  the  ray  making  a 
parabolic  trajectory  as  seen  in  Figure  4. a.  From  the  analytic  result  of  Budden,  we  know 
that  the  trajectory  in  this  ideal  linearly  varying  plasma  is  perfectly  parabolic. 

Table  1  compares  the  differences  between  the  trajectories  in  Figure  4.a  by 
examining  the  differences  between  the  y  coordinate  of  the  ray  trajectories  at 
corresponding  x  values.  The  analytic  result  of  Budden  and  the  numeric  results  of 
Haslegrove’s  equations  are  considered  to  be  the  correct  answer  because  they  are  the 
established  ray  tracing  methods. 

Table  1.  Comparison  of  Eikonal  Methods  in  a  Linearly  Varying  Plasma 


Category  _  Measurement 

Average  difference  between  y  coordinate  at  corresponding  x  coordinate  253.74  m 
for  the  SIM  and  Haslegrove’s  equations 

Average  difference  between  y  coordinate  at  corresponding  x  coordinate  254.28  m 
for  the  SIM  and  Budden  equation 

Average  difference  between  y  coordinate  at  corresponding  x  coordinate  0.55  m 
for  the  Haslegrove’s  equations  and  Budden’s  equation _ 
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—  Haslegrove  Differential  Equations 
—  Budden  Analytic  Solution 

b)  x(km) 


Figure  4.  Comparison  of  Eikonal  Methods  in  a  Linearly  Varying  Plasma 
a)  Trajectory  b)  Index  of  Refraction 

The  second  comparison  case  examined  the  trajectory  of  an  EM  wave  in  a  plasma 
exponentially  varying  in  density  in  one  dimension. 
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—  Haslegrove  Differential  Equations 
—  Budden  Analytic  Solution 

Figure  5.  Comparison  of  Methods  in  an  Exponentially  Varying  Plasma 

The  index  of  refraction  decreases  exponentially  in  the  y  direction,  which  results  in  a 
trajectory  very  similar  to  the  linearly  varying  plasma  case  except  that  the  trajectory  is  no 
longer  parabolic.  Figure  5  shows  that  the  trajectories  predicted  by  each  method  are  close 
enough  to  each  other  that  they  are  virtually  indistinguishable. 


Table  2.  Comparison  of  Eikonal  Methods  in  an  Exponentially  Varying  Plasma 


Category 

Measurement 

Average  difference  between  y  coordinate  at  corresponding  x  coordinate 
for  the  SEM  and  Haslegrove’ s  equations 

80.2  m 

Average  difference  between  y  coordinate  at  corresponding  x  coordinate 
for  the  SIM  and  Budden  equation 

79.7  m 

Average  difference  between  y  coordinate  at  corresponding  x  coordinate 
for  the  Haslegrove’ s  equations  and  Budden’ s  equation 

0.51  m 

A  symmetry  comparison  between  different  implementations  of  SIM  was 
performed  to  insure  that  the  SIM  was  properly  tracing  the  ray  path  of  the  EM  wave.  Due 
to  an  intuitive  understanding  of  Snell’s  law  and  the  analytic  results  of  Budden,  we  expect 
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the  path  of  an  EM  wave  refracted  by  a  plasma  with  perfectly  parallel  strata  to  refract 
through  the  atmosphere  symmetrically  (i.e.  if  the  trajectory  of  the  EM  wave  was  folded  in 
half,  the  halves  should  overlay  each  other  and  the  time  of  travel  should  be  the  same  for 
both  halves  of  the  trajectory).  The  following  graph  compares  the  symmetry  of  the 
calculated  trajectory  using  various  versions  of  SIM  for  an  EM  wave  traversing  a  linearly 
varying  plasma  with  parallel  strata. 

The  three  different  implementations  of  SIM  included: 

1)  Original  Algorithm  -  This  version  of  the  algorithm  represented  the  most  basic 
implementation  of  SIM.  It  is  simply  a  Euler  numeric  method  that  calculates  the  group 
velocity  and  radius  of  curvature  of  an  EM  wave  at  each  x,  y  coordinate  and  alters  the 
trajectory  of  the  EM  wave  according  to  the  magnitude  of  the  radius  of  curvature  at  a 
particular  point  in  the  plasma. 

2)  Predictor-Corrector  Methodology  -  uses  the  same  methodology  as  the  original 
algorithm  to  predict  the  next  point  of  the  trajectory  of  the  EM  wave.  The  P-C 
methodology  then  corrects  the  group  velocity  and  radius  of  curvature  by  averaging  their 
values  over  the  path  of  the  ray  and  uses  these  average  values  to  determine  the  next  point 
in  the  trajectory  of  the  ray. 

3)  Symmetric  Reflection  Algorithm  -  this  algorithm  insures  that  when  the  EM  wave  is 
reflected  in  the  plasma  that  the  reflection  is  symmetric  (i.e.  the  trajectory  of  the  EM  wave 

symmetric  about  the  vector  Vn ).  A  reflection  occurs  when  the  Z  component  of  the 

vector  resulting  from  the  cross  product  of  the  ray  direction,  s  ,  and  Vn  changes  sign. 
Symmetry  is  insured  by  transposing  the  position  and  angle  of  the  EM  wave  before  the 
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reflection  about  Vn .  This  is  done  so  that  the  angle  into  the  reflection  point  equals  the 
angle  exiting  the  reflection  point  resulting  in  a  symmetric  reflection. 


—  Predictor-Corrector 
--  Symmetric  Reflection 

Figure  6.  Symmetry  Comparison  of  the  Trajectory  of  an  EM  Wave 
in  a  Linearly  Varying  Plasma  a)  EM  Wave  Trajectory  b)  Comparison 
of  Relative  Difference  in  y  Coordinate  at  Corresponding  x  Values  over 
EM  Wave  Trajectory 

Figure  6  indicates  that  the  predictor-corrector  and  symmetric  reflection  algorithms 
make  the  trajectory  of  the  EM  wave  substantially  more  symmetric  which  increases  the 
accuracy  of  the  trajectory,  propagation  time,  path  length  of  the  EM  wave.  The  accuracy 
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of  the  calculation  of  those  quantities  will  be  critical  when  determining  the  attenuation  and 
refraction  of  the  EM  wave  over  very  short  distances  of  the  plasma. 

According  to  Sommerfeld,  as  a  plane  wave  traverses  a  medium  with  changing 
index  of  refraction  and  parallel  strata,  the  quantity  ft  sin  <9  should  remain  constant. 
Therefore,  a  check  to  insure  that  the  original  SIM  algorithm  was  maintaining  a  constant 
value  of  nsinO  was  performed  using  an  EM  wave  propagating  through  a  linearly 
varying  plasma.  The  following  figure  is  a  graph  of  ft  sin  6  over  the  trajectory  of  the  EM 
wave  in  Figure  4. 


Figure  7.  Check  for  a  Constant  Value  of  ft  sin  0  in  a  Medium 
with  Parallel  Strata 

It  should  be  noted  that  the  value  ft  sin  6  only  varies  in  value  by  0.001  over  the  entire 
trajectory  of  the  ray. 

This  chapter  presented  the  ray  tracing  methods  of  Haslegrove,  Budden,  and  SIM 
as  well  as  a  comparison  of  these  methods.  The  SIM  compared  well  to  the  established 
analytic  results  of  Budden  and  the  numeric  results  of  Haslegrove’ s  equations.  It  also  was 
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validated  by  the  symmetry  check  and  the  constant  n  sin  6  check.  SIM  also  is  capable  of 
ray  tracing  in  two  or  three  dimensions.  Therefore,  the  SIM  will  be  the  model  of 
refraction  used  to  determine  the  trajectory  of  an  EM  wave  as  it  traverses  the  electron 
beam  generated  plasma.  Since  we  have  developed  a  means  to  calculate  the  refraction  of 
an  EM  wave  in  a  plasma,  we  now  will  develop  the  equations  for  the  attenuation  of  an  EM 
wave  as  it  traverses  the  plasma. 


Attenuation  of  an  EM  Wave  in  a  Collisional  Plasma 

The  time-dependent  wave  equation  derived  from  Maxwell’s  equations  is 

k 2  d2E 


=  0 


c o 2  dt1 

where  the  plane  wave  solution  to  the  second-order  differential  equation  is  given  by 


(25) 


i(k-r) 


and 


where 


E  =  E0el(*' 


p  —  p  ,,-»(«+*) 

-c'o  —  •L'oe 


E0  =  electric  field  phasor 


E0  =  electric  field  amplitude 


The  wave  number  for  the  plane  wave  described  in  equation  (26)  is  given  by 


k  = 


net) 


(26) 


(27) 


(28) 


Since  the  index  of  refraction  is  complex,  the  wave  number  is  complex,  therefore 


E  =  EQe-{k‘~r)ei{K'7) 


(29) 
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where  the  term  e  (ki'r)  represents  the  amplitude  attenuation  of  the  beam  as  it  travels  a 
distance  r  through  the  medium.  The  attenuation  of  the  EM  wave  intensity  over  a  distance 

|r  I  is  expressed  by 


r, 


(30) 


where 

r7  =  intensity  attenuation 

Ia  =  intensity  of  the  attenuated  EM  wave 

I  =  intensity  of  the  unattenuated  EM  wave 
If  we  consider  a  plasma  whose  complex  wave  number  is  changing  as  a  function  of 
position  in  the  plasma  then  a  more  appropriate  attenuation  equation  is  given  by 
integrating  the  complex  index  of  refraction  over  the  path  of  the  ray  through  the  plasma. 
Integrating  over  the  path  of  the  ray  results  in  the  equation 

rl 

T,  =  exp(-2 -dr)  (31) 

r0 


which  describes  the  attenuation  of  an  EM  wave  over  its  entire  path  through  the  plasma. 
Equation  (31)  will  allow  us  to  determine  the  intensity  or  power  attenuation  of  an  EM 
wave  as  it  propagates  through  the  plasma.  In  the  next  section  equation  (31)  is  combined 
with  the  SIM  to  create  a  high  fidelity  model  of  the  refraction  and  attenuation  of  an  EM 
wave  in  a  plasma. 
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Implementation  of  the  EM  Wave  Propagation  Model  (EMWPM) 


To  estimate  the  attenuation  and  refraction  of  an  EM  wave  as  it  traverses  an 
artificially  generated  plasma,  a  wave  propagation  model  capable  of  computing  the 
trajectory  and  the  attenuation  of  an  EM  wave  at  any  frequency  in  a  two-dimensional, 
inhomogeneous  medium  was  required.  To  meet  this  requirement,  a  fortran  program  was 
written  which  combined  the  ray  tracing  method,  SIM,  the  amplitude  attenuation  model 
given  by  equation  (31),  and  a  two-dimensional  linear  interpolation  method  (described 
later  in  this  section). 

SIM  was  incorporated  into  EMWPM  by  using  a  Euler  predictor-corrector  method 
with  an  adaptive  step  size  (which  is  described  in  greater  detail  in  the  section,  Comparison 
of  Analytic  and  Numeric  Rav  Tracing  Results).  The  predictor  calculates  the  next  step 
using  the  ray  angle  and  group  velocity  at  the  current  point.  The  corrector  modifies  the 
ray  angle  and  group  velocity  by  performing  a  weighted  averaging  of  the  curvature  of  the 
ray  and  the  group  velocity  over  the  length  of  the  predictor  step,  hence  producing  a  more 
accurate  step.  If  the  distance  between  the  end  points  of  the  predictor  step  and  the 
corrector  step  are  larger  than  the  error  threshold  set  by  the  user  then  the  time  step  of  the 
calculation  is  halved  until  the  difference  between  the  end  points  is  within  the  error 
threshold. 

The  intensity  attenuation  of  the  EM  wave  is  based  on  equation  (31)  where  |Arj  |  is 
the  path  length  of  each  ray  segment  calculated  by  SIM  and  the  imaginary  part  of  the  wave 
number,  kt ,  is  given  by 

ki=ni-  (32) 

c 
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where  nj  is  calculated  from  Chapter  I,  equation  (12).  From  equation  (31),  the  following 
equation  was  developed  to  calculate  the  total  power  attenuation  of  the  EM  wave  as  it 
propagates  through  the  plasma 

r,  =  exp  -  2Yj \K  I Ai}  |  (33) 

The  input  to  the  EMWPM  consists  of  a  grid  of  densities  of  the  plasma  that  are  generated 
using  the  Electron  Beam  Simulation  (which  is  described  in  Chapter  IV).  To  reduce  the 
number  of  grid  points  required  to  accurately  sample  the  plasma  density,  a  two 
dimensional  linear  interpolation  method,  based  on  the  Taylor  series  expansion  of  a 
function  of  two  variables,  was  used  to  determine  the  plasma  density  between  grid  points. 
The  linear  interpolation  is  performed  using  the  following  equation 

5  =  s+^A  +  ^B  +  ^AB  (34) 

2  2  2 

where 


ml  +  m2  +  m3  +  m4 
s  — 

4 

(35) 

^  x-(x2+xl)/2 

{x2-xx)/2 

(36) 

B_  V  -  (V2  +  Vi)/2 

(y2-y  i)/2 

(37) 

m2+m2  ml  +  m4 
x~  2  2 

(38) 

A  mi+mi  m{  +  m2 

y  ~  2  2 

(39) 

A  ml+m3  m2  +  m4 

xy~  2  2 

(40) 
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m„  =  plasma  density  at  that  point 

x„  yj  =  coordinates  of  the  grid  end  points 

(Kiemele,  1997:8-15).  Figure  8  depicts  the  relationships  between  the  coordinates  (xj  ,yj) 
and  the  m„  variable.  The  terms  Ax  and  Ay  represent  the  change  in  the  plasma  density  in 

the  x  and  y  direction  and  the  term  represents  the  change  in  the  plasma  density  in  the 

diagonal  direction.  The  terms  A  and  B  are  the  coordinate  x  and  y,  respectively,  scaled  to 
a  value  between  -1  and  1.  Using  equations  (34)  -  (40),  the  plasma  density  can  be 
calculated  for  any  point  in  a  particular  cell.  As  the  EM  wave  propagates  through  the 
plasma,  EMWPM  determines  the  cell  in  the  grid  of  electron  densities,  which  is  required 
for  the  SIM  calculation.  EMWPM  then  performs  a  linear  interpolation  using  equations 
(34).(40)  to  determine  a  closer  approximation  of  the  plasma  density  along  the  path  of  the 
EM  Wave. 
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Figure  8.  Diagram  of  a  Cell  in  the  Plasma  Density  Table 
where  mi  is  the  Plasma  Density  at  the  Point  (xi,yi) 

The  SIM,  attenuation,  and  linear  interpolation  algorithms  described  in  this  section 
provide  the  core  of  the  EMWPM  program.  A  further  description  of  how  these  algorithms 
fit  together  is  provided  in  the  next  section. 

Description  of  Functions  and  Subroutines 

This  section  describes  the  main  algorithms  and  subroutines  as  well  as  the  logical 
flow  of  the  EMWPM  program.  For  a  top-level  flow  diagram  of  the  EMWPM  program 
see  Figures  9  and  10.  Table  3  contains  a  brief  description  of  the  functions  and 
subroutines  in  the  EMWPM  program.  The  SIM  method  described  in  the  previous 
subsections  is  incorporated  into  the  EMWPM  program  via  the  subroutine 
PlasmaRefractFunction.  The  linear  interpolation  method  described  in  the  previous 
section  is  implemented  in  the  NumberDensity  subroutine  and  is  essential  to  all 
calculations  of  the  plasma  frequency,  index  of  refraction,  group  velocity,  etc.  of  the  EM 
wave  in  the  plasma.  The  first  column  in  Table  3  provides  the  name  of  the  subroutine  or 
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function  in  the  EMWPM  program.  The  second  column  gives  a  brief  description  of  the 
function  a  particular  subroutine  followed  by  a  more  in-depth  description  of  the 
subroutine. 


Table  3.  Description  of  EMWPM  Subroutines  and  Functions 


Subroutine  Name 

Brief  Description 

MainProg 

Main  Program 

Obtains  input  from  the  user  from  the  Fortran  function 
NAMELIST  I/O,  checks  input  files  for  errors,  imports  the 
plasma  density  table,  calls  the  PlasmaRefractFunction,  and 
outputs  results  (see  Figure  9.).  Results  include  the  ray 
trajectory,  group  velocity,  index  of  refraction,  and  intensity 
attenuation  of  the  EM  wave. 

PlasmaRefractModule 

PlasmaRefractFunction 

Implements  numeric  solution  of  SEM  Differential  Eauations 

The  controlling  algorithm  for  the  EMWPM  program. 
PlasmaRefractFunction  both  controls  the  flow  of  the 

EMWPM  and  calls  all  the  functions  used  to  calculate  the 
refraction  and  attenuation  of  the  EM  wave  (see  Figure  10  for 
details). 

AirDensity 

Calculate  Nm 

Calculates  the  air  density  at  the  altitude  given  by  the  user. 
AirDensity  assumes  the  atmosphere  is  exponential  and  uses  a 
scaling  height  of  8180  m  (Al’pert,  1960:84) 

CollisionFreq 

Calculate  v 

Calculates  the  collision  frequency  between  thermal  electrons 
at  a  temperature  specified  by  the  user  and  neutral  air 
molecules  using  Chapter  I,  equation  (8). 

NumberDensity 

Calculate  Nc 

Determines  the  cell  that  a  coordinate  falls  in  by  dividing  the  x 
and  y  coordinate  of  interest  by  the  cell  width  and  height 
respectively.  This  provides  the  location  of  the  array  element 
in  the  three  dimensional  array  that  describes  the  plasma 
density  distribution.  Once  the  density  at  the  four  comers  of 
the  cell  is  established,  NumberDensity  function  linearly 
interpolates  using  equations  (34)  -  (40)  providing  an  estimate 
of  the  density  at  any  point  in  between  the  cell’s  comers 

PlasmaFreq 

Calculate  coF 

Calculates  the  plasma  frequency  using  Chapter  I,  equation  (1) 
based  on  the  number  density  provided  by  the  function 
NumberDensity. 
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GroupVelocity 

Calculate  v„ 

Calculates  the  group  velocity  of  the  plasma  using  Chapter  I, 
equation  (4). 

IndexofRefraction 

Calculate  nr 

Calculates  the  real  index  of  refraction  of  the  plasma  using 
Chapter  I,  equation  (11). 

WaveNumberlmag 

Calculate  kt 

Calculates  the  imaginary  wave  number  of  the  EM  wave  using 
equations  (8)  and  Chapter  I,  equation  (12). 

GradlndexofRefraction 

Calculate  Vn 

Calculates  the  gradient  of  the  index  of  refraction  of  the  EM 
wave  using  a  three-point  difference  formula  in  each  direction 
with  the  Ax  specified  by  the  user 

PropagationDirection 

Calculate  s 

Calculates  the  ray  propagation  unit  vector,  s  ,  in  Cartesian 
coordinates  using  the  equation  s  =  cos(#)x  +  sin((9)y . 

CurvatureVector 

Calculate  K 

Calculates  the  curvature  vector  describing  the  refraction  of 
the  EM  wave  as  it  propagates  through  the  plasma  using 
equation  (2) 

Magnitude 

Calculate  v| 

Calculates  the  magnitude  of  a  vector 

CrossProduct 

Calculate  v,  x  v. 

Determines  cross  product  of  two  arbitrary  vectors 

GetNew  Angle 

Calculate  Ad 

Computes  a  new  angle  of  propagation  for  the  ray  by  sampling 
the  plasma  density  at  the  end  points  and  the  mid  point  of  the 
predicted  path  of  the  ray.  A  new  group  velocity  and  radius  of 
curvature  is  calculated  for  each  of  these  points,  then  those 
values  are  averaged,  using  a  weighted  average  toward  the 
midpoint,  to  obtain  a  better  estimate  of  the  radius  of  curvature 
and  group  velocity  of  the  ray.  The  averaged  group  velocity 
and  radius  of  curvature  are  then  used  in  equation  (5)  to 
determine  the  change  in  angle  of  the  ray.  The  term  Viz  x  s  , 
determines  if  the  change  in  angle  is  added  or  subtracted  from 
the  initial  angle.  If  Vn  x  ?  is  positive  in  the  z  then  the  change 
in  angle  is  added  to  the  initial  angle,  and  if  Viz  x  s  is  negative 
the  change  in  angle  is  subtracted. 
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Figures  9  and  10  provide  an  overview  of  how  the  subroutines  described  in  Table  3 
work  together  in  the  EMWPM  program.  Figure  9  describes  the  data  flow  in  the 
EMWPM  program  for  a  time  varying  plasma.  The  logical  flow  of  the  subroutine, 
PlasmaRefractFunction,  is  presented  in  Figure  10  and  is  a  direct  result  of  the  modified 
Euler  method,  the  predictor-corrector  algorithm,  and  the  adaptive  time  step  that  are  used 
to  solve  the  SIM  differential  equations. 
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Figure  9.  MainProg  Flow  Diagram  from  the  EMWPM  Program 
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Figure  10.  PlasmaRefract  Function  Flow  Diagram  from  the  EMWPM  Program 
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EMWPM  Inputs  and  Outputs 


The  user  controls  the  simulation  EMWPM  via  a  Fortran  NAMELIST  file.  The 
file  contains  such  parameters  as  the  number  of  simulations  to  run,  maximum  dimensions 
of  the  plasma,  initial  starting  point  of  the  EM  wave,  number  of  rays,  spatial  extent  of  the 
modeled  wave,  average  electron  temperature  and  the  frequency  of  the  EM  wave.  The 
plasma  density  is  read  from  a  file  containing  the  x  and  y  position  as  well  as  the  plasma 
density  at  those  coordinates.  If  the  plasma  is  varying  in  time,  then  a  file  describing  the 
density  of  the  plasma  at  each  time  step  is  required  (This  is  acquired  from  Electron  Beam 
Simulation  described  in  Chapter  IV). 

EMWPM  outputs  a  single  file  for  each  simulation  run,  containing  data  on  the 
trajectory  of  the  ray,  propagation  time,  index  of  refraction,  Vn ,  group  velocity,  ray 
attenuation,  plasma  frequency,  and  complex  wave  number,  kt .  If  multiple  rays  are 

simulated,  the  file  is  divided  into  multiple  sections,  with  one  section  containing  the 
complete  history  of  one  of  the  simulated  rays.  If  multiple  simulations  are  run,  then 
multiple  output  files  are  produced  with  each  file  containing  the  ray  trajectory  histories  for 
a  certain  set  of  parameters.  If  the  plasma  density  varies  in  time  then  a  file  describing  the 
ray  path  for  each  time  step  is  produced  as  well. 

Capabilities  of  the  EMWPM  Program 

The  EMWPM  program  is  a  highly  flexible  beam  propagation  and  attenuation 
model  that  allows  the  user  to  simulate  multiple  rays  refracting  and  attenuating  in 
arbitrary,  inhomogeneous  plasma.  The  multiple  rays  may  be  used  to  represent  an  EM 
wave  of  finite  spatial  extent  refracting  and  attenuating  in  a  plasma.  The  EMWPM  is 
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capable  of  modeling  an  EM  wave  with  any  frequency  in  the  Radio  Frequency  (RF)  to 
optical  range,  but  EMWPM  assumes  that  the  wave  is  at  a  single  frequency  (i.e.  negligible 
bandwidth  for  purposes  of  attenuation  and  refraction).  EMWPM  is  also  capable  of 
modeling  an  inhomogeneous  plasma  that  varies  in  density  over  time. 

Now  that  we  have  a  program  for  determining  the  attenuation  and  refraction  of  the 
EM  wave  as  it  propagates  in  an  arbitrary  plasma,  we  must  develop  a  means  to  describe 
the  spatial  density  distribution  of  the  plasma.  Since  the  plasma  will  be  generated  through 
electron  impact  with  air  molecules,  the  next  section  will  describe  electron  impact  theory, 
which  includes  both  elastic  scattering  and  electron  impact  ionization  of  neutral 
molecules. 
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III.  Electron  Impact  Cross  Section  Theory 

To  develop  a  code  to  calculate  the  spatial  extent  and  density  distribution  of  the 
electron  beam  generated  plasma,  the  collision  cross  sections  between  an  electron  and  the 
molecules  in  the  air  must  be  obtained.  Due  to  limited  experimental  data  on  electron 
impact  ionization  cross  sections  at  electron  energies  between  2  keV  and  5  MeV,  we  must 
use  theoretical  models  to  obtain  the  ionization  cross  sections  for  our  plasma  generation 
model.  Currently  there  are  several  theories  that  have  been  developed  to  describe  the 
cross  section  of  an  electron  colliding  with  a  neutral  molecule.  The  elastic  scattering  cross 
sections  discussed  in  this  section  deal  with  the  electron  scattering  due  to  the  coulomb 
field  of  the  nucleus.  The  ionization  cross  sections  deal  with  a  free  electron  colliding  with 
an  atomic  or  molecular  electron.  Some  of  these  models  (such  as  Mott’s  ionization  cross 
section)  not  only  describe  the  probability  that  the  incident  electron  will  ionize  the 
molecule,  but  also  the  energy  lost  and  the  angle  scattered  by  the  incident  electron  as  well 
as  the  energy  and  direction  of  the  ejected  electron.  This  information  can  be  used  in  a 
simulation  that  models  the  trajectories,  energy  loss,  and  number  of  ionizations  that  occur 
as  a  beam  of  electrons  propagates  through  the  air.  From  the  results  of  such  a  simulation, 
we  will  be  able  to  determine  the  density  of  the  plasma  generated  by  a  relativistic  beam  of 
electrons  ionizing  the  air. 

Background  Theory  of  Scattering 

Rutherford  scattering  provides  a  classical  view  of  how  an  electron  is  scattered  by 
another  charged  particle  due  to  the  interaction  of  the  electron  with  the  coulomb  field  of 
the  nucleus.  Consider  a  fast  electron  passing  near  a  nucleus  of  charge  Ze  and  mass  M. 
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The  perpendicular  distance  between  the  electron  velocity  and  the  nucleus  is  referred  to  as 
the  impact  parameter,  b.  From  Figure  1 1,  we  see  that  as  the  electron  passes  by  the 
nucleus  it  is  deflected  by  an  angle  d  due  to  the  coulomb  attraction  between  the  electron 
and  the  nucleus. 


According  to  Fermi, 

“The  cross  section  for  scattering  of  the  incident  particle  at  an  angle  6  in 
the  range  dd  is  defined  to  be  the  total  area  perpendicular  to  the  initial 
path  of  the  particle  such  that  if  the  particle  passes  through  this  area  it  is 
deflected  by  an  angle  0  in  dd  ” 


(Orear,  1949:35).  The  area  perpendicular  to  the  initial  velocity  of  the  electron  that  will 
scatter  into  the  angle  6  in  the  range  dd  given  by 

dff(d)  =  27ib(d)db(d)  =  2nb(d)^^-dd  (1) 

dd 


Classical  mechanics  gives  the  following  formula  for  the  relationship  between  the  angle  of 
deflection  and  the  ratio  of  the  potential  and  kinetic  energy  for  two  particles  interacting  via 
a  coulomb  force  (for  details  of  the  derivation  see  (Evans,  1955:843)) 


d  Ze 2 

tan— = - t- 

2  myb 


(2) 


the  impact  parameter,  b,  being  specified  as: 
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I//1N  Ze2  0 
b(0)  = - -cot— 


(3) 


my 


We  can  calculate 


db(6) 

dd 


for  purposes  of  substituting  into  equation  (1) 


db(0)  Ze 2  2  0 

— —  = - -esc  — 

d0  2my2  2 


(4) 


giving  us  a  single  differential  cross  section  (SDCS)  described  by 

da{0)  =  - 


7tZ2e*  l2  0  2  0  rn 

-b  cot— esc  —d0 


(5) 


(mev2)2  2  2 

To  determine  the  expression  for  a  particle  scattered  into  a  solid  angle  Q.  in  the  range  of 
dQ ,  we  can  easily  change  equation  (5)  using  an  expression  for  a  differential  solid  angle 

(6) 


0  0 

d£l  =  2n  sin  0d0  =  An  sin— cos— d0 

2  2 


resulting  in  the  equation 

dcr{0)  = 


ZV 


1 


■da 


(7) 


4 (mev2)2  sin4  0/2 

Hence,  we  have  developed  the  nonrelativistic,  SDCS  for  the  elastic  scattering  of 
an  electron  through  the  solid  angle  dQ  in  the  laboratory  frame  of  reference.  To  convert 
to  a  relativistic,  differential  cross  section  the  following  relationships  can  be  used  for  the 
mass  and  velocity  of  the  electron 

v  =  j3c  (8) 


m. 


m  = 


Vw?5 


(9) 


which  results  in  the  relativistic  differential  cross  section 
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(10) 


dt 7(0)  = 


z2zV  (i— fi2 


2  \  2 


4(m  c  ) 


ytf4  )  sin4  0/2 


<£2 


(Evans  1955:593).  This  equation  results  in  a  scattering  angle  distribution  for  an  electron 
scattering  off  a  bare  nucleus  of  an  atom.  By  integrating  (10),  we  can  obtain  the  total 
elastic  scattering  cross  section  of  an  electron  with  a  bare  nucleus.  However,  the  equation 
has  a  singularity  at  a  scattering  angle  of  zero  and  from  experiment  we  know  that  the 
probability  of  the  electron  scattering  into  the  angle  6  -  0  is  not  infinite.  Therefore,  a 
common  method  of  circumventing  the  flaws  in  this  classical  approach  is  to  use  a  lower 
limit  of  integration  other  than  0  (Lawson,  1988:257).  The  value  of  the  lower  limit  is 
discussed  in  the  next  section,  which  presents  a  quantum  mechanical  approach  to 
determining  the  elastic  scattering  angle  distribution  of  the  electron. 


Mott’s  Elastic  Scattering  Cross  Section 

Due  to  the  Heisenberg  uncertainty  principle,  the  classical  theory  of  scattering  is 
limited  to  the  domain  where 


-2*  »i 

131  ft 


(ID 


(Evans,  1955:593).  Thus  Rutherford’s  scattering  cross  section  is  limited  to  slow  electrons 
colliding  with  a  nucleus  containing  many  protons.  Since  the  electrons  that  we  are 
interested  in  will  have  energies  on  the  order  of  106  eV  resulting  in  f5  =  1  and  a  nuclei, 
with  a  Z  of  7  or  8,  therefore,  a  quantum  mechanical  treatment  of  the  elastic  scattering 
cross  section  will  be  required. 

Using  the  relativistic  Dirac  theory  of  the  electron  and  the  First  Bom  Approximation,  Mott 
obtained  the  full  form  of  the  relativistic  differential  cross  section  of  particles  scattering 
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under  the  influence  of  a  coulomb  field.  He  then  developed  an  approximate  form  yielding 
a  differential  cross  section  very  similar  to  the  Rutherford  scattering  cross  section: 


dt 7(0)  = 


z2ZV  1  -p 


2 


4 (mec  )  l  fi 


0  0  0 

1  -  ft1  sin 2  —  +  n/3a(\  -  sin  — )  sin  — 
2  2  2 


sin4  0j 2 


dQ 


(12) 


where  the  new  term  represents  the  effect  of  electron  spin  and  indistinguishability  on  the 
scattering  of  the  electron  (Mott,  1965:235).  Integrating  equation  (12)  gives  the  total 
collision  cross  section,  however  equation  (12)  predicts  that  a  singularity  occurs  at  6  =  0 . 
From  experiment  we  know  that  the  probability  of  the  electron  scattering  into  the  angle 
0  =  0  is  not  infinite.  Therefore,  a  common  method  of  circumventing  this  problem  is  to 
use  a  lower  limit  of  integration  other  than  0  (Lawson,  1988:257).  The  minimum  angle  of 
scattering  for  an  electron  scattering  off  a  nucleus  will  occur  when  the  electron  impact 
parameter  is  approximately  the  same  as  the  atomic  electron  screening  radius.  According 
to  Lawson,  0^  can  be  calculated  by 


Q  X  _  Xar2Z^  a 


rmn  A 


Py 


\meJ 


(13) 


(Lawson,  1988:257),  where 

X  =  de  Broglie  wavelength  of  the  incident  particle 


r 

rs  =  — £—r  =  electron  screening  radius 
a2Z^ 

m0=  rest  mass  of  the  electron 

The  maximum  angle  of  scattering  for  an  electron  scattering  off  a  nucleus  will  occur  when 
the  electron  undergoes  a  head  on  collision  with  the  nucleus  and  therefore  is  Pi,  except  in 
the  ultra-relativistic  cases  (Lawson,  1988:257).  Therefore,  the  total  elastic  cross  section 
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for  an  electron  scattering  off  a  nucleus  after  integrating  equation  (12)  is  described  by  the 
equation 

4^e4Z2 
2  c4/3*me2y2 

Ionization  Cross  Section 

Currently  there  are  several  theories  that  have  been  developed  to  describe  the  cross 
section  of  ionization  for  an  electron  impacting  a  neutral  molecule.  All  of  these  models 
can  be  used  to  obtain  the  total  ionization  cross  section  of  an  electron  impacting  a  neutral 
molecule.  However,  some  of  these  models  (such  as  Mott’s  ionization  cross  section) 
provide  additional  information  such  as  the  angular  scattering  distribution  of  the  incident 
electron  as  well  as  the  angle  and  energy  distribution  of  an  ejected  electron.  Using  these 
distributions,  we  can  develop  a  simulation  that  models  the  trajectories,  energy  loss,  and 
number  of  ejected  electron  produced  as  an  electron  travels  through  the  air.  From  those 
results,  we  will  be  able  to  determine  the  density  of  our  plasma  due  to  the  electron  beam 
firing  into  the  air. 

Bethe’s  Relativistic  Ionization  Cross  Section 

Bethe  performed  a  detailed  quantum  mechanical  calculation  using  the  First  Bom 
Approximation  to  determine  the  average  energy  lost  by  a  fast  particle  when  colliding 
with  an  electron  bound  to  a  nucleus.  His  perturbation  calculation  starts  with  the  coulomb 
potential  between  the  bound  electron  and  the  stationary  nuclear  charge,  Ze.  He  then 
added  two  perturbation  terms  that  represent  the  potential  energy  between  the  incident 


cxtiZJS  (esc ^2—  1)  +  i (esc2  ^^--1)  + 

a(7tZ(3  +  —)  In  j sin  ^s. . 
a  2 


(14) 
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particle  and  the  nucleus  and  the  bound  electron.  His  solution  was  extended  to  atoms 
containing  Z  electrons  by  replacing  the  standard  coulomb  potential  of  a  bare  nucleus  with 
a  field  due  to  a  bare  nucleus  plus  the  field  due  to  the  (Z  -  1)  atomic  electrons  (Evans 
1955:579).  This  results  in  a  non-relativistic  ionization  cross  section  that  is  proportional 
to  l/<22,  where  Q  is  the  energy  of  the  ejected  electron,  and  an  energy  loss  per  unit  path 
length  described  by 


dT_ 

ds 


4  7Z2e4 


iVZ  In 


f  2  m0V2^ 


V 


(15) 


where 

T  =  kinetic  energy  of  the  incident  electron 
V  =  velocity  of  the  incident  particle 
mo  =  rest  mass  of  the  incident  particle 
z  =  charge  of  the  incident  particle 

I  is  the  geometric  mean  of  all  the  ionization  and  excitation  potentials  of  the  atom 
involved  in  the  collision.  I  is  defined  as 

=  (16) 

^  n,l 

where  fn ,  is  the  sum  of  the  oscillator  strengths  for  all  optical  transitions  of  the  electron  in 

the  n,  l  orbital  and  is  on  the  order  of  unity  for  most  atoms.  An  l  is  the  mean  excitation 

energy  of  the  n,  l  orbital  and  its  value  is  fairly  close  to  the  ionization  energy  of  the  orbital 
electrons  for  the  outer-shell  electrons.  However,  theoretical  values  of  7  for  atoms  other 
than  hydrogen  are  hardly  ever  used  because  they  are  usually  incorrect.  Therefore,  the 
experimentally  determined  values  of  7  are  commonly  used  (Evans  1955:579). 
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Bethe  modified  equation  (21)  to  account  for  the  Lorentz  contraction  of  the  electric 
field  of  a  relativistic  incident  particle.  This  general  result  is  applicable  to  soft  collisions 


where  the  energy  of  the  ejected  electron  is  between  (2mm  and  H,  where  Qmi„  is  the 
minimum  ejected  electron  energy  and  H  is  the  maximum  ejected  electron  energy 
considered.  This  modification  results  in  the  equation 


dT,  2*V  ,„,i  2m0V2H 


_ s__  _ 

ds  m0V 2 


NZ  In 


/2(1-/?2) 


~P2 


(17) 


(Evans,  1955:582).  Evans  states  that  equation  (23)  can  be  extended  to  all  electron  impact 
ionization  collisions  by  making  H  =  <2max  =  772.  This  only  results  in  an  error  of  a  few 
percent  from  a  more  exact  expression  for  energy  lost  per  unit  path  length.  Hence, 
equation  (17)  becomes 


dT  lire  , 

—  = - -NZ  ln| 

ds  m0V 


4  f  mnV2T 


I2(l-j32) 


-P2 


(18) 


A  cross  section  can  then  be  obtained  by 


a  = 


dT  1  d  (Ionization) 


N- 


dT 


ds  N 


ds 


(19) 


d  (Ionization) 


which  results  in 


a  = 


2 7te 


m0V 4 


-In 


avg 


m0Vzr 

I2(l-J32) 


-P2 


(20) 


where 


Iavg  =  average  ionization  energy  of  an  atom  or  molecule 
The  following  graph  displays  the  total  ionization  cross  section  of  molecular  nitrogen 
obtained  from  equation  (20)  versus  the  incident  electron  energy 
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Figure  12.  Total  Ionization  Cross  Section  from  Bethe’s  Equation 
Mott’s  Ionization  Cross  Section 


Mott  also  obtained  an  electron  impact  ionization  cross  section  for  a  fast  electron 
impacting  an  electron  bound  by  a  generalized  coulomb  potential.  Using  the  First  Bom 
Approximation,  the  differential  cross  section  of  a  free  electron  colliding  with  a  bound 
electron  is  obtained.  The  differential  cross  section  of  that  collision  is  given  by  the 
expression 


ImK(d)dKdco  = 


4 n2m2  kn 
~hi  k 


-|JjV  exp\i(km/cnl  -  kn0)  ■  R\f/*Ky/mdrdR  dK da)  (21) 


where 


V  =  coulomb  interaction  energy  between  the  incident  and  atomic 


electrons  e2  / 


R 


k  =  incident  electron  wave  number 
k  =  wave  number  of  the  continuous  spectrum  state 
m  =  initial  primary  quantum  number  of  the  bound  electron 
kmK  =  the  wave  number  of  the  incident  electron  after  ionization 
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y/m  =  bound  electron  wave  function 

y/K  =  electron  wave  function  in  the  continuous  spectrum 

r  -rx-r2  =  distance  between  the  atomic  and  incident  electron 

R  =  ^-(rl  +r2  )  =  coordinate  of  the  center  of  mass  of  the  electrons 

Mott  solved  this  equation  using  a  wave  function,  developed  by  Sommerfeld,  for  an 
electron  in  the  continuous  spectrum,  K ,  moving  in  a  direction  corresponding  with  polar 
angles  ■  The  Sommerfeld  wave  equation  is  given  by 

y/*K  =y/K(r,n-®)  =  (  In)  *  K  exp  {iter  + — n^)|  u~me~u  J0  \l-yj  im^]du  (22) 

^  o 


where 

£  =  r(  1  +  cos  0) 

cos  0  =  cos  Q  cos  ^  +  sin0  sin  ^f(cos  <j>  cos  yr  +  sin  <f>  sin  y/) 


n  =  Z  /  Ka0 


(Mott,  1965:489).  Mott  obtained  an  ionization  TDCS  for  a  fast  incident  electron 
colliding  with  an  atomic  or  molecular  electron  with  a  particular  binding  energy  by 
assuming  that  the  effects  of  interference  between  the  ejected  and  incident  electrons  is 
small  and  hence  negligible.  Using  this  form  of  the  wave  equation  (21),  Mott  obtained  the 


differential  cross  section 

2 S/J6ke  ks  exp[-(2// !ke) tan -1  {2/zfce  7(//2  +Ak2  -fcj} 


IK  ( G)dad(odK 


m02Nc2  k  (1-e  2nfl,k‘)(ju2  +  Ak2  +ke2  -2AkkecosSY 


x 


(A k  -  ke  cos  S )2  +  ju2  cos2  8 
(ju2  +Ak2  -k2)2  +4ju2k2 


(23) 


dcrdcodK 


where 
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//  —  /  a0 


A k  = 


A 

ksks  -kk 


S  =  the  angle  between  ksks  -  kk  and  the  direction  described  by 

(%,y/)  =  electron  ejection  angles  relative  to  the  incident  electron  direction 

dc r  =  differential  solid  angle  in  which  the  electron  is  ejected 

(0,0)  =  electron  scattering  angle  relative  to  the  incident  electron  direction 

dco  =  differential  solid  angle  in  which  the  incident  electron  is  scattered 

ke  =  ejected  electron  wave  number 

ks  =  scattered  electron  wave  number 

ao  =  Bohr  radius 

Figure  13  summarizes  the  relationships  between  the  wave  vectors  in  this  equation. 
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The  differential  cross  section,  Equation  (23),  obtains  its  maximum  value  when  S  =  0 , 
which  corresponds  to  momentum  being  conserved  among  the  electrons.  The  maximum 
of  equation  (23)  is  given  by  the  conservation  of  momentum  condition 

k2  +  ke2  -  2kke  cos x  =  K2  (24) 


(Mott,  1965:490).  Mott’s  equation  assumes  that  the  velocity  of  the  incident  electron  is 
much  larger  than  the  velocity  of  the  bound  electron,  thus  making  the  kinetic  energy  of  the 
bound  electron  negligible.  Therefore,  the  energy  conservation  condition  can  be  written  as 

T  =  W  +  TS  +  B  (25) 


where 


and 


T  =  incident  electron  energy 
W  =  ejected  electron  energy 
Ts  =  scattered  electron  energy 
B  =  binding  energy  of  the  atomic  electron 

n2k2 

T  =  (26) 

2  me 


Therefore, 

k2=ke2+ks2+k02  (27) 

From  equation  (23),  the  angular  distribution  of  the  scattered  electron  can  be  found  by 
integrating  over  all  angles  of  ejection.  This  can  be  accomplished  by  obtaining  an 

expression  for  cos  <5  from  the  scalar  product  of  ke  and  A k  which  results  in  the 


expression; 
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(28) 


k-Ak  = 


A  k 


cos< 


Solving  equation  (28)  for  cost)  in  spherical  coordinates  yields 

cos  S  =  cos  6  cos  x  +  cos  </)  sin  6  cos  y/  sin  x  +  sin  V sin  X  sin  0  sin  </>  -  cos  x  (29) 

The  resulting  analytic  equation  after  integration  is 

,  ,/n  2 wju6ke  ks  exp[-(2/// ke )  tan-1  {2//fce  /(//2  +Ak2  -*,)}„ 

IK{G)d(odK=  2;— — - (1_e-2^/^} - x 


a^AV 
A k2+±(ke2+ju2) 


(30) 


| «4  +2/^2(M2  +A:£2)  +  (A^2  ~ke2)2} 


■dcodtc 


(Mott,  1965:490).  Using  equation  (30)  and  (23),  the  scattering  angle  distribution  of  the 
incident  electron  and  the  angular  and  ejection  energy  distribution  of  the  ejected  electron 
can  be  determined  for  the  collision  of  a  free  electron  with  an  electron  in  a  generalized 
coulomb  field.  It  should  be  noted  that  equations  (23)  and  (30)  indicate  that  the  ejected 
and  scattered  angle  distributions  are  isotropic  in  y/  and  <j)  respectively.  Figures  14  and  15 
give  an  example  of  an  ejected  and  scattered  electron  angular  distribution  for  an  electron 
in  the  3(7g  orbital  of  molecular  nitrogen. 

As  the  incident  electron  energy  increases,  small  scattering  angles  become  very 
favored,  such  that  the  scattering  angle  distribution  peaks  near  0  =  0.  This  indicates  that 
at  high  electron  energies,  forward  scattering  of  the  incident  electrons  is  highly  favored. 

As  a  result  of  higher  incident  electron  energies,  the  ejection  angle  distribution  in  Figure 
14.b  spikes  at  the  angle  corresponding  to  the  conservation  of  momentum. 
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Figure  14.  Mott’s  Angular  Distributions  for  Molecular  Nitrogen 
for  an  Incident  Electron  of  100  eV  and  Ejected  Electron  of  20  eV 
a)  Scattering  Angle  b)  Ejection  Angle 


The  ejected  electron  energy  distribution  is  calculated  by  integrating  equation  (30) 
over  all  ejection  angles.  This  integration  must  be  done  numerically  and  results  in  an 
ejection  energy  distribution  seen  in  Figure  15. 
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Figure  15.  Mott’s  Ejection  Energy  Distribution 
for  Molecular  Nitrogen  at  an  Incident  Electron 
Energy  of  80  eV 

It  should  be  noted  that  upon  integrating  Mott’s  TDCS  (equation  (13))  three  times 
that  the  total  cross  section  for  a  single  electron  in  a  particular  orbital  can  be  calculated. 
To  obtain  the  total  ionization  cross  section  for  the  entire  molecule,  the  cross  section  for 
each  electron  in  each  molecular  orbital  must  be  summed  together.  This  method  of 
calculating  the  total  cross  section  of  ionization  is  very  computationally  intense 
(considering  that  the  simulation  will  need  to  calculate  this  quantity  several  million  times 
to  determine  if  an  electron  experiences  an  elastic  or  ionization  collision),  therefore  an 
analytic  equation  is  needed  to  obtain  the  total  ionization  cross  section  of  a  molecule. 

Binary  Encounter  Bethe  Ionization  Cross  Section 

Kim  and  Rudd  of  NIST  developed  a  SDCS,  which  they  called  the  Binary- 
Encounter-Dipole  model  (Kim,  2000:052710-1).  It  uses  a  modified  form  of  the  Mott 
cross  section  for  the  collision  of  two  electrons  in  the  presence  of  a  generalized  coulomb 
potential 


57 


4 m2R2N 


1 


1 


1 


(1 W  +  B)2  (W  +  B)(T-W)  (T-WY  J 


(31) 


where 


R  =  the  Rydberg  energy 
N  =  orbital  occupation  number 

This  equation  does  not  take  into  account  dipole  interactions,  which  are  soft 
collisions  that  result  in  a  small  transfer  of  momentum  to  bound  electrons.  Therefore,  the 
modified  Mott  cross  section  is  combined  with  the  Bethe  cross  section  for  soft  collisions 
to  obtain  the  BED  equation.  Kim  and  Rudd  required  that  the  combined  Mott  and  Bethe 
formula  satisfy  asymptotic  forms  for  both  the  ionization  and  stopping  cross  sections  of 
Mott  and  Bethe.  This  requirement  succeeded  in  eliminating  empirical  parameters  that 
had  been  used  in  previous  attempts  to  combine  the  two  equations.  The  SDCS  form  of  the 
BED  is  given  as 


B(t  +  u+ 1) 
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where 


S  =  4mQ2N(R/B)2 

(33) 

t-T  IB 

(34) 

u-U  IB 

(35) 

w-W  I B 

(36) 

N,=  jf^V 

(37) 
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—  =  differential  dipole  oscillator  strength 
dw 


U  =  orbital  energy  of  the  bound  electron 

To  obtain  the  relativistic  BED  equation,  Kim  and  Rudd  convert  the  nonrelativistic 
electron  velocities  and  energies  into  their  relativistic  counterparts: 

P  =1l  /? 2  =  1 - l—  t'=Tlmc 2  (38) 

c  1  +  t'2 

A=-  A2  =1-7 -W  i>'=  Blmc1  (39) 

c  1  +  b 

B  =^JL  B2=  1 - l—  u'-  U /me2  (40) 

Pu  c  Pu  1  +  u'2 


where 


v,  =  the  relativistic  speed  of  an  electron  with  energy  T 
vb  =  the  relativistic  speed  of  an  electron  with  energy  B 
vb  =  the  relativistic  speed  of  an  electron  with  energy  U 
Substituting  relationships  (38)  to  (40)  into  equation  (32)  we  obtain  the  following 
equation: 
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(41) 


By  integrating  equation  (41)  over  w  from  0  to  (t  - 1)/  2 ,  we  obtain  an  expression  for  the 
total  ionization  cross  section.  The  limit  (t  —  1)  /  2  is  the  maximum  amount  of  energy  that 
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is  given  to  the  ejected  electron  due  to  the  assumption  that  the  electron  leaving  the 
collision  volume  with  the  most  energy  is  the  incident  electron.  Since  the  differential 

dipole  oscillator  strength  term,  — ,  is  not  always  known  for  a  molecule  it  is 

dw 


approximated  by  Kim  and  Rudd  using  a  simple  function  that  simulates  the  shape  of 


dw 


in  the  ionization  of  the  hydrogen  atom.  In  the  case  when  no  data  is  available  for  the 
dipole  constant,  Kim  and  Rudd  set  the  dipole  constant  of  the  molecule  equal  to  one  (Kim, 
2000:0527103).  This  approximation  results  in  the  relativistic  Binary-Encounter-Bethe 
model,  given  by  the  expression 
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(42) 


The  relativistic  form  of  this  equation  is  required  for  electron  energies  greater  than  20 
keV,  however  this  formula  reduces  to  the  non-relativistic  BEB  equation  in  the  low  energy 
limit  (Kim  2000:0527101).  This  equation  describes  the  ionization  cross  section  for 
atomic  or  molecular  electrons  in  a  given  orbital,  and  thus  allows  us  to  be  specific  about 
which  molecular  electron  that  the  incident  electron  ionizes.  To  obtain  the  total  ionization 
cross  section  for  a  molecule,  the  cross  sections  for  each  shell  are  added  together. 

The  reason  that  three  ionization  cross  sections  are  reviewed  in  this  section  is 
because  they  all  have  their  strengths  and  weaknesses  for  being  used  in  the  electron  beam 
model.  Bethe’s  cross  section  requires  minimal  information  and  computation  (average 
ionization  energy,  I,  from  equation  (16),  and  Z)  to  provide  a  total  ionization  cross  section 
for  a  molecule.  Mott’s  equation  (23),  however,  provides  the  ability  to  calculate  the 
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scattering  angle  distribution  of  the  incident  electron  as  well  as  the  ejected  angle  and 
energy  distributions  of  the  liberated  electron.  Mott’s  equation  (23)  allows  us  to  develop  a 
sophisticated  model  of  the  electrons  propagating  through  the  air.  Mott’s  equation  needs 
to  be  numerically  integrated  three  times  in  order  to  obtain  a  total  ionization  cross  section, 
which  is  overly  computationally  intense  for  a  Monte  Carlo  simulation.  Therefore,  the 
RBEB  model  is  a  good  complement  to  Mott’s  equation  because  it  provides  an  analytic 
solution  for  the  total  ionization  cross  section  of  an  orbital  shell  and  has  been  extensively 
validated  with  experimental  data.  Therefore,  a  combination  of  the  two  models  was  used 
in  the  Monte  Carlo  simulation  that  is  described  in  Chapter  IV.  In  the  simulation,  Mott’s 
elastic  cross  section  model  and  the  RBEB  model  calculate  the  total  elastic  and  ionization 
cross  sections  respectively,  which  enables  the  Monte  Carlo  simulation  to  determine  if  the 
electron  experiences  an  elastic  or  ionization  collision.  If  the  electron  undergoes  an 
ionization  collision,  then  Mott’s  ionization  equation  provides  the  scattering  and  ejected 
angle  distribution  as  well  as  the  ejection  energy  distribution,  which  determines  the  energy 
and  direction  of  the  scattered  and  ejected  electrons  after  the  collision  with  the  air 
molecule. 

Comparison  of  Cross  Section  Results 

A  comparison  of  Mott’s  differential  elastic  cross  section,  equation  (12),  with  data 
from  NIST’s  elastic  scattering  data  base  was  performed  at  incident  electron  energies  of 
50, 10,000,  and  20,000  eV  for  atomic  nitrogen.  Since  no  data  has  been  found  comparing 
Mott’s  equation  (12)  to  experimental  data  and  NIST  has  extensively  validated  their  model 
with  experimental  results,  it  is  appropriate  that  we  compare  the  more  simplistic  model  of 
Mott  to  the  NIST  model  developed  by  Jablonski,  et.  al. 
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•  NIST  Data  60 

_  Mott  Equation 

Figure  16.  Comparison  of  Mott  and  NIST  differential  cross  sections 
a)  50  eV  b)  10,000  eV  and  c)  20,000  eV 

From  Figure  16. a,  we  see  a  large  difference  between  the  Mott  and  NIST  differential  cross 
sections.  This  is  due  to  the  first  Bom  approximation  no  longer  being  valid  when  a  slow 
electron  (50  eV)  interacts  with  a  small  nucleus.  However,  Figure  ll.b  and  ll.c  show  that 
at  greater  incident  electron  energy,  the  first  Bom  approximation  becomes  valid.  Hence, 
the  NIST  and  Mott  differential  cross  sections  are  very  similar  at  high  incident  electron 
energies.  Since,  the  electrons  we  are  concerned  with  are  in  the  greater  than  10  keV  range 
(electrons  with  energies  less  than  10  keV  do  not  travel  very  far  through  the  air  (<1  cm)), 
therefore,  Mott’s  differential  cross  section  equation  is  adequate  for  describing  the 
scattering  angle  distribution  of  the  electrons  at  these  higher  energies. 

A  comparison  between  NIST  ionization  cross  section  data  and  the  relativistic 
Bethe,  Mott,  and  RBEB  total  ionization  cross  section  of  nitrogen  was  also  performed. 

For  the  Mott  and  RBEB  ionization  cross  sections  only  the  first  four  orbital  shells  were 
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used  in  the  computation.  This  is  because  the  inner  shells  do  not  contribute  significantly 
to  the  ionization  cross  section  and  also  RBEB  and  Mott’s  equations  do  not  accurately 
model  the  ionization  of  inner  shell  electrons  (Kim,  2000:052710-10).  The  molecular 
orbital  constants  for  nitrogen  are  given  in  Table  4.  To  calculate  Mott’s  ionization  cross 
section,  an  effective  Z  must  be  obtained  for  the  variable// .  Since  an  effective  Z  is  due  to 
the  shielded  coulomb  field  that  is  acting  on  the  electrons  in  their  molecular  orbital,  the 
measured  binding  energies  were  used  to  calculate  an  effective  Z  using  the  equation 

Z  =  (43) 

B  =  binding/ionization  energy  of  the  electron 


Table  4.  Molecular  Orbital  Constants  of  Nitrogen 


Molecular  Binding  Energy  (eV)  Average  Electron  Dipole 

Orbital  Kinetic  Occupation  Constant 

Energy  Number 

(eV)  ‘ 


lCTg 

427.41 

601.38 

2 

1 

lau 

427.30 

602.40 

2 

1 

2ag 

41.72 

71.13 

2 

1 

2<7u 

21.00 

63.18 

2 

1 

lTTu 

17.07 

44.30 

4 

1 

3crg 

15.58 

54.91 

2 

1 

Remarks:  Data  was  obtained  from  (Kim,  2000) 

Binding  energy  from  experimental  vertical  ionization  energy 


For  Bethe’s  equation  an  average  ionization  energy  of  36.6  eV  (Evans,  1955:659)  and  I 
value  of  86  eV  (Evans,  1955:583)  were  used  in  the  calculation.  It  was  unclear  what  the 
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appropriate  value  for  Z  should  be  for  molecular  nitrogen,  therefore  the  two  extreme 
values  of  7  and  14  were  used  to  determine  the  possible  range  of  Bethe’s  total  ionization 
cross  section.  This  determination  of  the  extreme  values  for  Bethe’s  total  ionization  cross 
section  is  shown  in  Figure  17.a. 


Mott  Electron  Energy  (eV) 

—  NIST 
•  RBEB 


Figure  17.  Comparison  of  Bethe,  RBEB,  Mott,  and  NIST  Ionization  Cross  Sections 
a)  Comparison  of  RBEB  and  Bethe’s  Total  Ionization  Cross  Sections  b) 
Comparison  of  NIST,  Mott’s,  and  RBEB’s  Total  Ionization  Cross  Sections 
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RBEB,  Mott’s  and  NIST  cross  sections  appear  to  be  in  very  good  agreement  at  lower 
energies  with  a  very  small  difference  between  their  values.  Whereas  Bethe’s  cross  section 
marginally  agrees  with  the  NIST/RBEB  model  depending  on  the  value  of  Z  used  in  the 
calculation  and  the  energy  range  over  which  the  cross  section  is  compared.  At  higher 
incident  electron  energies,  Bethe’s  cross  section  provides  a  fairly  close  estimate  to  the 
total  ionization  cross  section. 

In  this  chapter,  theories  on  electron  elastic  and  ionization  collisions  were 
reviewed.  The  mathematical  collision  cross  section  models  resulting  from  the  theories  of 
Mott,  Bethe,  and  Kim  were  discussed  to  gain  insight  into  their  capabilities  and 
limitations.  The  resulting  cross  sections  from  each  model  were  compared  to  cross  section 
values  that  have  been  established  by  NIST.  The  next  chapter  will  show  how  these  cross 
section  models  are  utilized  in  a  Monte  Carlo  model  that  will  determine  the  plasma  density 
distribution  created  by  an  electron  beam. 
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IV.  Electron  Beam  Propagation 


Simple  Electron  Beam  Propagation  Model  (SEBPM) 

The  primary  emphasis  of  the  SEBPM  is  to  develop  a  lower  and  upper  bound  to 
the  longitudinal  extent  of  the  plasma  generated  by  an  electron  beam.  The  SEBPM  is 
limited  to  a  scenario  where  the  electron  travels  through  the  air,  but  experiences  no 
angular  deflection  in  a  collision.  Two  limiting,  energy  partitions  associated  with 
ionization  are  considered:  the  ejected  electron  is  stationary  and  the  incident  electron 
energy  is  reduced  by  the  ionization  energy  or  the  scattered  and  ejected  electrons  share  the 
incident  electron  energy,  reduced  by  the  ionization  energy,  equally.  The  electrons  are 
then  propagated  through  the  air  until  all  the  electrons  have  energies  less  than  the  average 
ionization  energy  of  the  air  molecules.  These  ultimate  electrons  are  considered  thermal 
electrons.  The  limited  scenarios  described  above  were  addressed  to  provide  a  bound  on 
the  spatial  distribution  of  the  plasma  and  relate  these  results  to  those  derived  from  the 
more  complex  Monte  Carlo  simulation. 

Axial  Density  Profile  with  No  Angular  Scattering 

The  scenarios  described  previously  were  named  1)  Ejected  Electrons  Receive  No 
Energy  (RNE)  scenario  and  2)  Ejected  and  Scattered  Electrons  Equally  Share  Energy 
(ESE)  scenario.  In  an  ionizing  collision,  the  incident  electron  identifies  the  electron 
before  the  collision,  the  scattered  electron  is  the  incident  electron  after  the  collision,  and 
the  ejected  electron  is  the  electron  that  is  liberated  from  the  target  molecule  or  atom.  In 
the  RNE  scenario,  the  incident  electron  has  a  relativistic  velocity  and  collides  with  a 
bound  electron,  which  results  in  an  ionization  event  where  one  of  the  electrons  (either 
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incident  or  ejected)  leaves  the  collision  volume  (defined  by  the  collision  parameter  b) 
with  the  incident  energy  minus  the  binding  energy  of  the  molecular  electron.  The  ejected 
electron  leaves  the  collision  volume  as  a  thermal  electron.  It  is  assumed  that  all 
ionization  events  in  the  RNE  scenario  result  in  only  the  loss  of  the  average  ionization 
energy  from  the  incident  electron.  The  RNE  scenario  establishes  the  maximum 
longitudinal  extent,  because  the  incident  electron  is  neither  scattered  nor  is  any  energy 
given  to  the  ejected  electron. 

The  ESE  scenario  is  similar,  but  it  assumes  that  the  ejected  and  scattered  electrons 
both  emerge  from  the  collision  volume  with  an  equal  amount  of  energy.  This  type  of 
ionization  event  is  termed  by  Evans  as  a  “Hard  Collision”  and  they  are  infrequent  and 
therefore,  contribute  very  little  to  the  most  probable  energy  loss  of  an  electron.  The  ESE 
scenario  is  the  other  bound  of  the  possible  length  of  the  electron  beam  generated  plasma, 
because  the  electrons  reduce  in  energy  by  half  after  every  collision  they  do  not  travel 
very  far  through  the  air. 

Model  Theory 

The  following  differential  equation  describes  the  decay  of  a  beam  of  particles  due 

to  single  collisions  experienced  by  the  particles  as  the  beam  propagates  in  the  x  direction 

dp_  _  P  /j\ 

dx~  A 

The  solution  to  this  equation  represents  an  exponential  loss  in  the  intensity  of  the  beam 
over  a  distance  x 

P  ~  P initial*'  ^  (2) 

where 
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p  -  number  density  of  the  beam  of  particles  (#/cm  ) 

However,  if  a  single  collision  does  not  stop  the  particle,  but  rather  lowers  the  energy  of 
the  particle,  we  may  use  a  similar  differential  equation  to  describe  the  loss  of  the 
electrons  from  that  energy  state. 

dPp  _  A) 
dx  A0 

where 

p0  =  density  of  the  beam  of  particles  at  the  initial  energy 

A0  =  mean  free  path  (MFP)  of  particles  at  the  initial  energy  of  the  beam 

If  we  make  the  assumption  that  the  particle  losses  the  average  ionization  energy  of  the 
gas  every  time  it  undergoes  an  ionizing  collision,  then  we  can  discretize  the  energy  states 
such  that  we  have  energy  states  0, 1,  2  ...  n.  Where  the  nth  energy  state  is  defined  as  the 
energy  of  an  electron  after  it  has  experienced  n  ionizing  collisions.  Therefore,  the  energy 
of  the  nth  energy  state,  En ,  is  given  by 

E„  =  E0 

where 

E0  =  initial  energy  of  the  electron 
A  E  =  average  ionization  energy  of  the  gas 

From  this  assumption  we  can  write  the  equations  for  the  change  in  density  of  particles  in 
the  lower  energy  states  as 

d£i  _  Po  _  Pi  ^4) 

dx 
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(5) 


_  Pi  Pi 

dx  ^2 


to 


dPn  =  Pn-1  Pn 
dx  Xn_x  Xn 


where 


X.  =■ 


1 


iW7(E„) 


(7) 


N  =  number  density  of  the  air  molecules 
cr(En)=  ionization  cross  section  as  a  function  of  the  incident  electron  energy 


Here,  <J(En)  is  calculated  from  the  Bethe,  Mott,  or  BEB  ionization  cross  section 

equation.  However,  this  results  in  an  enormous  number  of  equations  for  the  RNE 
scenario.  As  an  example,  a  beam  of  1  MeV  electrons  would  result  in  27,322  energy 
states  above  thermal  energy  and  hence  would  result  in  27,  322  differential  equations  to 
solve! 

Therefore,  it  was  necessary  to  develop  a  general  analytic  solution  for  the  system 
of  first  order  differential  equations  (4)-(6).  For  the  RNE  scenario  the  general  solution  to 
the  system  of  equations  described  by  equations  (4)  through  (6)  is 


p  =f  (-1  YN&'Xs 


X 

At 


k=0 


m 


(8) 


;= o 

j*k 


N  =  initial  number  of  electrons  in  the  electron  beam 
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For  the  ESE  scenario,  equations  (4)  to  (6)  are  modified  by  a  factor  of  2”  because  two 
electrons  of  the  same  energy  are  created  in  each  collision  therefore 

d-Pn  _  _  P_n 

dx  An_x  An 


which  results  in  a  general  solution 


X 


=  ^(r\y2-NX'\e± 

r* n  /  a 


Jfe=0 


(10) 


j= 0 
j*k 

For  the  RNE  scenario  a  stationary,  thermal  electron  is  created  in  every  ionization 
collision,  as  a  result  the  number  of  electrons  created  is  determined  by  summing  the 
electron  densities  in  all  pn  states.  The  result  of  the  summation  will  yield  the  over  all 


axial  electron  density  profile  of  the  electron  beam.  For  the  ESE  scenario,  each  electron 
emerged  from  the  collision  volume  with  equal  energy,  therefore  the  initial  energy  of  the 
electron  is  roughly  halved  after  every  collision  resulting  in  a  minimal  number  of  electron 
energy  states  (for  a  1  MeV  electron  there  would  be  14  energy  states).  Therefore,  the 
general  solution  (10)  can  be  used  to  calculate  the  electron  number  profile  for  the  ESE 
scenario  fairly  easily. 

However,  there  is  a  numerical  problem  associated  with  the  RNE  scenario  and  the 
use  of  equation  (8).  When  and  Ak  are  very  close  in  value  there  is  a  loss  in  numerical 

precision,  which  increases  with  each  state  due  to  the  product  term  in  the  denominator  of 
equation  (8).  This  loss  in  precision  eventually  leads  to  completely  inaccurate  results. 
Therefore,  a  numerical  method  that  replicates  the  physics  of  equation  (8)  was  developed 
to  estimate  the  electron  number  profile  of  the  electron  beam. 
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Numeric  Model  Theory 


The  numeric  method  for  the  SEBPM  is  implemented  by  treating  each  energy  state 
as  a  separate  beam  of  electrons  passing  through  the  medium.  In  the  RNE  scenario  an 
ionization  event  results  in  the  same  loss  of  energy  each  time,  therefore  the  number  of 
electrons  in  each  beam  decays  exponentially  with  distance.  Equation  (2)  represents  the 
number  of  electrons  that  do  not  experience  a  collision  in  a  distance  x.  Therefore: 

P  =  P„m(  I-'-*)  <11) 

describes  the  number  of  electrons  that  experience  a  collision  in  a  distance  x.  The  numeric 
method  or  cascading  method  works  by  demoting  the  electrons  that  underwent  a  collision 
to  a  lower  energy  state  each  time  a  collision  occurs.  Those  electrons  that  do  not  undergo 
a  collision  stay  in  the  same  energy  state.  This  process  is  repeated  until  each  electron 
experiences  enough  ionizing  collisions  that  its  energy  is  less  than  the  average  ionization 
energy  of  the  gas.  Figure  18  illustrates  how  the  electrons  decay  from  one  energy  state  to 
another  as  the  electron  beam  travels  through  a  distance,  Ax .  The  electrons  that  decay 
from  a  higher  energy  state  provide  a  source  to  the  electron  beam  in  the  next  lower  energy 
state  (see  Figure  18  for  a  pictorial  description  of  this  process).  This  numeric 
approximation  is  very  similar  to  the  analytic  form,  but  does  not  suffer  from  the  loss  in 
numeric  precision.  The  cascading  method,  however,  is  very  computational  intense  when 
the  electrons  have  a  large  initial  energies  because  of  the  large  number  of  energy  states 
that  must  be  tracked  and  the  large  number  of  Ax  intervals.  For  the  cascading  method  of 
SEBPM,  the  value  of  Ax  must  be  approximately  the  same  or  less  than  the  mean  free  path 
of  the  electrons  at  that  particular  energy  state.  If  the  Ax  value  was  set  much  larger  than 
the  MFP  of  the  electrons,  then  all  the  electrons  decayed  to  the  lower  energy  state  in  one 
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Ax  step.  This  resulted  in  an  overestimation  of  the  longitudinal  extent  of  the  electron 
profile,  because  the  electrons  may  have  experienced  multiple  collisions  with  in  the  large 
Ax  step,  but  the  model  only  allowed  one  ionization  collision  to  occur.  If  Ax  was  set  to 
small,  then  the  numerical  cascade  model  took  a  very  long  time  to  execute  because  it  had 
to  perform  many  steps  to  calculate  the  longitudinal  propagation  distance  of  the  electron 
beam.  Therefore,  the  Ax  step  had  to  be  reasonably  close  to  the  smallest  electron  MFP 
value  in  the  simulation  to  obtain  a  reasonable  estimate  of  the  longitudinal  extent  of  the 
plasma. 


Eo 

Ei 


E2 


En 


Figure  18.  Diagram  of  Energy  States  in  the 
SEBPM  Numerical  Method 

Results  of  Simple  Electron  Beam  Propagation  Model 

For  the  SEBPM  additional  simplifying  assumptions  were  introduced  to  make 

calculations  easier.  First,  the  air  is  assumed  to  be  100%  molecular  nitrogen  and  second, 

an  average  energy  of  36.6  eV  is  assumed  to  be  lost  per  ionization  event  (Evans, 
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1955:698).  Bethe’s  total  ionization  cross  section,  Chapter  III,  equation  (20),  was  used  in 
all  calculations  with  an  average  ionization  energy  of  36.6  eV,  an  I  of  86  eV,  and  a  Z  of 
14. 


Results  of  the  RNE  Scenario 

In  Table  5,  the  results  of  the  Cascading  Method  for  the  RNE  scenario  at  various 
initial  electron  energies  are  shown.  The  power  of  the  electron  beam  is  fixed  at  500  kW 
and  the  electron  beam  is  operated  for  1  second.  Hence,  as  the  initial  electron  energy 
increases,  the  initial  number  of  electrons  decreases  proportionally.  The  Ax;  for  the  model 
was  set  at  0.01  cm  which  is  close  in  magnitude  to  the  minimum  mean  free  path,  at  the 
lowest  energy  state  in  the  model.  Table  5  shows  the  results  of  running  the  model  with 
initial  electron  energies  of  1  MeV,  2  MeV  and  5  MeV  at  a  constant  power  of  500  kW. 
The  relationship  between  the  power  of  the  electron  beam,  P,  and  the  number  of  electrons 
in  the  electron  beam,  Ne ,  is 


P 


Me 

t 


N„  = 


EL 

E0 


where 


E0  =  energy  of  the  electrons  in  the  electron  beam 
Table  5.  Results  of  Cascade  Model  for  a  RNE  Scenario 


Electron  Energy 

Starting  #  of 
Electrons 

Total  #  of  Electrons 
after  Ionization 

Range  of  Electron 
Profile 

1  MeV 

3.121xl01!i 

8.529xl022 

1639  cm 

2  MeV 

1.561xl018 

8.528xl022 

3579  cm 

5  MeV 

6.242xlOn 

8.530xl022 

9100  cm 
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The  resulting  longitudinal  density  profile  for  an  electron  beam  at  1  and  2  MeV  is 
shown  in  Figure  19.  The  peak  values  in  the  longitudinal  density  profiles  corresponds  to 
the  peak  value  of  the  ionization  cross  section  for  the  lower  energy  electrons  (see  Figure 
17.a). 


Figure  19.  Electron  Longitudinal  Density  Profile  from  the  RNE  Scenario 
for  an  Electron  Beam  of  a)  1  MeV  b)  2  MeV 

Results  of  the  ESE  Scenario 

The  following  table  shows  the  results  from  equation  (10)  at  various  initial  electron 
energies.  The  power  of  the  electron  beam  is  fixed  at  500  kW  and  the  pulse  length  of 
electron  beam  is  1  second.  For  this  scenario  Ax  is  set  at  0.01  cm,  which  is  close  to  the 
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magnitude  of  the  mean  free  path  of  the  lowest  energy  state.  An  assumption  was  made 
that  the  ionization  cross  section  for  ESE  scenario  was  the  same  as  the  RNE  scenario  for 
an  equivalent  energy  of  the  incident  electron.  This  assumption  is  not  true,  but  it  did 
present  the  extreme  limiting  case  for  the  distance  that  the  electron  beam  would  travel. 

The  final  electron  distribution  is  not  found  the  same  way  as  in  the  RNE  scenario,  because 
the  ejected  electrons  have  the  same  energy  as  the  scattered  electrons.  Therefore,  the 
electron  distribution  profile  is  the  electron  number  of  the  nth  or  last  energy  state.  From 
equation  (6),  we  obtain  the  relationship  between  the  pn  and  pn_x  which  is 


Pn 


=  [  —  ^P-dx' 

J  ; 

0  1 


(12) 


The  term  An_x  has  a  large  range  of  possible  values  depending  on  the  final  energy 

state  in  the  calculation.  This  is  due  to  the  cross  section  peaking  as  the  incident  electron 
energy  nears  the  binding  energy  of  the  molecule  or  atom.  After  peaking  the  ionization 
cross  section  decreases  rapidly  in  value  (see  Figure  17).  The  value  of  pn  is  very 
dependent  on  the  value  of  An_x ,  and  since  An_x  has  such  a  large  range  of  possibilities;  we 
must  insure  that  we  obtain  a  reasonable  value  for  An_x .  Therefore,  two  methods  were 
implemented  of  determining  An  ,  for  the  ESE  scenario.  In  the  first  method,  the  mean  free 
paths  of  the  last  five  energy  states  (An_x,An_2  ,An_3 . ..)  are  averaged  over  an  energy  bin. 
Since,  the  difference  in  energy  between  the  En.j,  En. 2,  and  En.3 ,  corresponding  to 
A  ,  ,A  ,,A  ,  for  the  ESE  case  is  not  constant,  therefore  the  energy  bin  is  defined  as  the 

n— 1 7  n—l  7  n—j 

energy  half  way  between  the  next  highest  and  lowest  energy  states. 

A£  (EnZEn-x)  +  (En-l  ~  En-l)  (13) 
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For  the  second  method,  no  averaging  is  done  for  the  Xn_x .  The  results  of  the 

electron  distribution  for  both  methods  and  the  length  of  the  electron  beam  profile  are 
given  in  Table  6  and  the  resulting  electron  number  density  profiles  are  shown  in  Figures 
20, 21,  and  22. 


Table  6.  Results  of  ESE  Scenario  for  Initial  Energies  of  1  MeV,  2  MeV,  and  5  MeV 


Electron 

Energy 

1  MeV 

2  MeV 
5  MeV 


Starting  #  of 
Electrons 

3.121xlOre 

1.561xl018 


6.242x10* 


#  Electrons 
Ionized  with 
X  averaged 
5.11  lx  1022 
5.1136xl022 


4.0909x10 


Electrons 
Ionized  with 
X  not  averaged 
6.184xl022 
6.187xl022 
4.267xl022 


Range  of 
Electron  Profile 

0.7  cm 
0.8  cm 
1.2  cm 


(a)  Averaged  MFP  (b)  Non- Averaged  MFP 

Determining  if  the  electrons  conserve  energy  provides  a  simple  check  to  the  validity  of 
the  models.  The  total  ending  energy  is  calculated  by  multiplying  the  total  number  of 
electrons  made  by  the  average  ionization  energy  for  the  RNE  scenario.  Comparing  the 
total  ending  energy  with  the  total  starting  energy  of  the  electron  beam  we  obtain  Table  7. 


Table  7.  Test  for  Energy  Conservation  for  the  RNE  scenario 


Electron  Energy 

Total  Starting  Energy  (J) 

Total  Ending  Energy  (J) 

1  MeV 

5x10s  J 

500,086  J 

2  MeV 

5xl05 J 

500,040  J 

5  MeV 

5x10s J 

500,134  J 

Table  7  indicates  that  the  model  conserves  energy  almost  perfectly.  For  the  ESE  scenario 
the  total  ending  energy  is  calculated  by  multiplying  the  total  number  of  electrons  made  by 
the  energy  in  the  last  energy  state.  Comparing  the  total  ending  energy  with  the  total 
starting  energy  of  the  electron  beam  we  obtain  Table  8,  which  indicates  that  the  ESE 
model  conserves  energy  almost  perfectly  as  well. 
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Table  8.  Test  for  Energy  Conservation  for  the  ESE  scenario 


Electron  Energy 

Total  Starting  Energy 

(J) 

Last  Energy  State 

Total  Ending  Energy 

(J) 

1  MeV 

”  5x10*7“ 

61  eV 

4.99x10s  J 

2  MeV 

5x10s J 

61  eV 

4.99xl05  J 

5  MeV 

5x10s  J 

76  eV 

4.98xl05  J 

The  SEBPM  bounded  the  electron  beam  longitudinal  extent  for  various  initial 
electron  energies  (See  Table  9.).  Because  the  power  attenuation  of  the  wave  is  equal  to 
e~2k,r ,  we  need  the  plasma  to  have  both  large  spatial  extent  and  density  (because  k, 

increases  with  plasma  density,  Chapter  II,  equation  (4».  The  maximum  electron  beam 
longitudinal  extent  is  approximately  proportional  to  the  initial  electron  energies  of  the 
electron  beam,  if  we  assume  that  the  transverse  extent  of  the  plasma  stays  the  same  with 
increasing  electron  energy,  then  the  density  of  the  plasma  will  reduce  linearly  with 
increasing  electron  energy.  If  we  assume  that  the  transverse  extent  of  the  electron  beam 
will  also  increase  linearly  with  the  initial  electron  energy  (which  is  shown  to  be  the  case 
in  Chapter  V),  then  the  decrease  in  plasma  density  will  be  proportional  to  1/T3 ,  where  T 
is  the  initial  energy  of  the  electrons  in  the  electron  beam.  Therefore,  the  initial  electron 
energies  of  the  electron  beam  will  most  likely  be  at  1  MeV  or  below. 


Table  9.  SEBPM  Minimum  and  Maximum  Beam  Length  Results 


Electron  Energy 

Min  Electron  Beam 
Penetration 

Max  Electron 

Beam  Penetration 

Total  Power 

1  MeV 

0.7  cm 

1639  cm 

500  kW 

2  MeV 

0.8  cm 

3579  cm 

500  kW 

5  MeV 

1.2  cm 

9100  cm 

500  kW 
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The  SEBPM  has  bounded  the  length  and  the  number  density  of  the  plasma,  but 
we  still  have  no  estimate  on  the  transverse  extent  of  the  electron  beam  due  to  electron 
scattering.  To  obtain  an  estimate  of  the  width  and  distribution  of  the  electron  beam 
generated  plasma,  a  Monte  Carlo  method  will  be  used.  The  Monte  Carlo  method 
described  in  the  next  section  of  this  chapter  should  give  a  more  accurate  result  as  to  the 
length,  width,  and  density  distribution  of  the  plasma  created  by  the  electron  beam. 

Monte  Carlo  Method 

The  purpose  of  this  section  is  to  describe  the  Monte  Carlo  method  and  show  how 
it  can  be  utilized  to  develop  a  more  detailed  model  than  was  used  in  the  previous  section. 
The  last  subsections  of  this  section  are  devoted  to  describing  the  program  that  was 
developed  to  simulate  the  electron  beam  propagating  through  the  air  using  the  Monte 
Carlo  method. 

The  Monte  Carlo  Method  models  random  processes  such  as  particle  diffusion  and 
transport  by  tracing  the  histories  of  sample  particles  as  they  travel  through  the  medium. 
All  collision  events  between  particles  are  assigned  a  probability  of  occurring  and  a 
pseudo-random  number  is  picked  to  determine  which  event  took  place.  The  results  of 
the  randomly  selected  collision  events  result  in  a  unique  trajectory  for  each  particle. 

After  calculating  the  histories  of  many  particles,  we  can  treat  the  results  of  the  simulation 
as  a  statistical  sample  of  how  all  the  particles  in  a  system  may  behave.  Therefore,  a 
correlation  must  be  established  between  statistical  and  physical  results  of  the  Monte 
Carlo  Simulation. 
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Monte  Carlo  Techniques 


The  purpose  of  this  section  is  to  provide  a  brief  overview  of  the  basics  of  the 
Monte  Carlo  method.  For  a  more  detailed  discussion  of  the  Monte  Carlo  Method  the 
reader  should  consult  the  Computational  Methods  of  Neutron  Transport  (Lewis 
1984:296-356).  In  the  case  of  the  relativistic  electron  beam  traversing  the  air,  the  focus 
will  be  exclusively  on  the  interactions  between  electrons  and  air  molecules.  In  Chapter 
TTT,  several  relationships  between  the  energy  of  an  electron  and  its  collision  cross  section 
with  a  molecule  were  developed.  If  we  were  to  introduce  an  electron  with  mass  me, 
velocity  ve,  charge  e  into  a  gas  of  molecules  of  number  density  Nm  with  velocity  «  ve,  in 
which  the  collision  cross  section  of  the  electron  with  the  molecules  is  <7,we  could 
describe  the  frequency  of  the  electron’s  collisions  with  the  air  molecules  by  the 
relationship 

v  =  Nma\ve\  (13) 


If  we  consider  many  electrons  in  the  gas  following  the  same  path,  we  can  describe  the 
rate  at  which  the  electrons  experience  a  collision  by  the  differential  equation 


(14) 


which  results  in  a  solution  that  represents  the  number  of  electrons  that  experience  a 
collision 


P  =  p0  exp[-J  vdt] 


(15) 


where 


p  =  density  or  intensity  of  the  electrons  in  the  beam  (e'/cm3  or  g/  cm3) 
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The  collision  frequency  is  constant  with  respect  to  t  since  there  are  no  forces  acting  on 
the  electrons,  therefore 

p  =  pQ  exp[-«f]  (16) 

The  exponential  term  represents  the  probability  density  function,  f(t),  that  the  electron 
will  experience  a  collision  between  time,  t  and  t  +  dt.  The  cumulative  probability 
distribution  function  is  defined  by 

F(t)  =  P{t'<  t)  (17) 


and  is  the  probability  that  the  random  variable,  t' ,  is  less  than  or  equal  to  t.  The 
relationship  between  the  probability  density  function  and  the  cumulative  probability 
distribution,  F(t),  is 


m= 


dF(t) 

dt 


(18) 


Therefore,  the  number  of  electrons  that  collide  in  the  time  interval  0  <t  <T  can  be 
represented  by  the  equation 

p  =  p0F(T)  =  p0  (1  -  P)  =  p0  (1  -  exp [-vt])  (19) 


Using  equation  (18),  we  can  also  introduce  the  rules  for  transformation  of  random 
variables.  We  begin  by  letting 

x  =  x(t)  (20) 

where  t  is  a  random  variable.  If  g(x)dx  is  the  probability  that  x  is  between  x  and  x  +  dx, 
and  f(t)dt  is  the  probability  that  t  is  between  t  and  dt  and  if  these  probabilities  are  equal 
then  we  can  write  the  relationship 

\g(x)dx\  =  \f(t)dt\  (21) 


Since  the  probability  distribution  functions  must  be  positive,  we  obtain  the  expression 
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g(x)  =  f(t)^- 
ax 


(22) 


If  we  consider  jc  being  the  particular  function,  x  =  Fit)  of  the  random  variable  t,  where 
F(t)  is  the  cumulative  probability  distribution.  As  a  result  equation  (22)  becomes 

Using  equation  (18)  to  determine  the  derivative  in  equation  (23),  we  obtain  for  the 
transformation,  x  =  F(t) ,  the  equation 

g(F(t))  =  1  (24) 

where 

g{F{t))dF{t)  =  dF{t)  (25) 

0  <  F{t)  <  1  (26) 

Thus  the  probability  of  the  random  variable,  Fit) ,  taking  on  a  value  between  F  and  F  + 
dF  is  equal  to  dF.  This  shows  that  F  is  uniformly  distributed  between  zero  and  one. 

A  Pseudo-Random  Number  (PRN)  from  a  computer  is  also  evenly  distributed  between  0 
and  1  and  is  unbiased.  By  setting  the  cumulative  probability  distribution  equal  to  the 
PRN  generated  by  a  computer 

F(T)  =  C  (27) 

where  £  is  a  random  variant  from  a  PRN  generator,  we  are  able  to  obtain  an  unbiased 
distribution  of  the  F(T)  values  from  the  computer’s  PRN  generator.  However,  we  are 
interested  in  the  distribution  of  t  values,  therefore  we  must  invert  the  cumulative 
distribution  function: 

t  =  F~\C)  (28) 
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For  the  simple  case  of  no  external  forces  on  the  electron  this  is  not  complicated,  but  for 
cases  when  there  is  an  external  force  on  the  electron  this  process  maybe  difficult.  The 
associated  cumulative  distribution  function  for  the  probability  density  function  described 
by  equation  (16)  is 


F(t)  =  £  =  l-exp[-Mf,] 


(29) 


The  result  of  the  inversion  for  equation  (29)  is 


(,=-±111(1-0  P°) 

V 

From  equation  (30),  we  obtain  a  distribution  of  the  time  between  collisions  for  the 
electrons.  With  the  time  distribution  and  the  electron’s  energy,  the  history  of  the  electron 
can  be  followed.  In  two  dimensions,  the  trajectory  of  the  electron  between  collisions  will 
be  simply 


X  =  |ve|*COS0 

(31) 

y  =  |ve|fsin0 

(32) 

where  0  will  change  after  each  collision.  If  we  start  with  a  simplifying  assumption  that 
after  a  collision  all  electrons  are  scattered  isotropically,  then  we  will  have  all  the 
equations  we  need  to  write  a  simple  Monte  Carlo  simulation  (Lewis,  1990:299-303). 

If  a  joint  probability  density  function  is  separable  such  that 

f(x,y)  =  Mx)f2(y)  (33) 


then  jc  and  y  are  said  to  be  independent.  For  an  isotropic  angular  distribution,  the  joint 
angular  probability  density  function  is  separable  as  well. 


sindd0d<f> 
4 n 


_  sin  6ddd</) 
4 n 


(34) 


where  C  is  a  constant.  Since 
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•sin  dddd<j) 
An 


=  1 


(35) 


then  the  following  equation  is  true 


/«?,*)  =  C=1 


(36) 


then 


sin  6d 0d<p  _(dAY dB" 
An  1  2  I  2n . 


(37) 


where  A  =  cos  6  and  B  =  </>.  Hence,  the  cumulative  distribution  simply  becomes 

F(A, B)  =  F(A)F(B)  Jl~A^  B 


(38) 


V  2  A  2  n) 

And  an  isotropic  scattering  angle  is  determined  from  the  equations  using  the  random 
variates 

F(A)  =  £ 

F(B)  =  Ci 


which  results  in  the  ith  scattering  angle  being 


0i  =  cos  (2£j  - 1) 


(39) 


<f>,  =  2 n£2  (40) 

Using  this  angular  distribution,  we  can  then  follow  the  electron  as  it  travels  through  the 
air.  For  the  two-dimensional  case  these  equations  will  enable  us  to  trace  the  electron 
between  collisions  using  the  following  equations 

xi  =  xi- 1  +  K  K  cos  0,-i  (41) 


=  A-i  +  A  t,  sin  0._x 


(42) 


where 
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xi ,  yt  =  the  position  of  the  electron  at  the  ith  collision 
=  time  between  the  i  and  i  - 1  Collisions 

and  0,_t  =  the  scattering  angle  in  laboratory  frame  after  the  i  - 1  collision 

(Ramos,  1990:46-47).  For  the  case  of  a  relativistic  electron  scattering  after  an  ionization 
or  elastic  collision,  an  isotropic  approximation  to  the  scattering  angle  distribution  is  not 
valid.  Therefore,  the  scattering  angle  distribution  developed  by  Mott,  Chapter  III, 
equation  (30),  should  be  used  for  the  Monte  Carlo  simulation. 

When  there  are  multiple  types  of  collisions  being  modeled  by  the  Monte  Carlo 
simulation  the  cross  section  can  be  interpreted  as  the  probability  of  having  a  collision  of  a 
given  type.  The  total  cross  section,  awml ,  is  equal  to  the  sum  of  all  possible  collision 

cross  sections. 

=2>,  («> 


where 


<ji  =  cross  section  for  a  collision  of  type  i 

Because  the  collisions  are  mutually  exclusive  events,  the  probability  that  a  particular  type 
of  collision  will  occur  is 


(44) 

^ total 

If  we  select  a  random  number,  g ,  that  is  uniformly  distributed  between  zero  and  one,  we 
can  determine  which  type  of  collision  occurs.  Hence,  we  can  determine  if  an  electron 
incident  on  a  molecule  experiences  an  elastic  or  ionization  collision  with  either  an 
oxygen  or  nitrogen  molecule  based  on  their  cross  sections  and  relative  concentrations. 
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For  further  discussions  on  the  Monte  Carlo  Method  the  reader  should  consult 
Computational  Methods  of  Neutron  Transport  (Lewis,  1984:296-303). 

Using  the  equations  presented  in  this  section  and  the  theoretical  cross  sections 
developed  in  Chapter  III,  we  can  develop  a  computer  program  that  tracks  the  trajectory  of 
the  electron  until  it  no  longer  has  enough  energy  to  ionize  another  molecule.  We  also 
note  that  each  ionization  event  results  in  another  electron  whose  trajectory  must  be 
tracked  as  well.  The  coordinates  of  the  thermalized  electrons,  electrons  that  have  energy 
less  than  the  minimum  ionization  energy,  are  recorded.  The  thermalized  electrons 
coordinates  are  then  used  to  determine  the  density  distribution  of  the  plasma  using 
descriptive  statistics,  which  are  presented  in  the  next  section. 

Descriptive  Statistics 

Because  there  will  be  a  large  number  of  ionization  events  per  high  energy 
electron,  even  if  we  start  the  Monte  Carlo  simulation  with  a  modest  number  of  initial 
electrons  we  could  end  up  with  millions  of  ejected  electrons.  If  a  1  Mev  electron  loses  all 
its  energy  in  ionization  collisions,  the  electron  will  have  liberated  27,322  molecular 
electrons  (given  an  average  energy  loss  of  36.6  eV  per  ionization).  If  we  start  the  Monte 
Carlo  simulation  with  1000  electrons  then  we  will  end  up  with  over  27  million  electrons. 
Tracking  and  performing  statistics  on  several  million  individual  electrons  requires  an 
exceptional  amount  of  computer  resources.  Therefore,  thermalized  electron  coordinates 
were  binned  in  a  two-dimensional  grid  and  group  statistics  were  used  to  determine  the 
distribution  characteristics  -of  the  plasma.  Group  statistics  are  used  when  the  data  has 
been  put  into  bins  and  individual  data  points  are  not  known.  The  following  equation 
gives  the  sample  mean,  x  ,  of  a  distribution  of  binned  data 
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(45) 


Z/A 

x  - 

n 

where 

k  =  number  of  bins  in  the  frequency  distribution 
Mi  =  midpoint  of  the  i,h  bin 
fi  =  frequency  of  the  ith  bin 
n  =  total  sample  size 

and  the  sample  variance  is  estimated  by  the  equation 

s2=4=! - - -  (46) 

ft  - 1 

(Kiemele,  1997:75-76).  The  Monte  Carlo  simulation  estimates  the  mean  and  standard 
deviation  of  the  electron  distribution,  transverse  to  the  axis  of  the  beam,  using  these 
equations. 


Electron  Beam  Simulation  (EBS)  (Monte  Carlo)  Description 

The  EBS  uses  a  Monte  Carlo  methodology  to  determine  the  electron  distribution 
resulting  from  a  relativistic  beam  of  electrons  propagating  through  the  atmosphere.  The 
simulation  was  written  in  Fortran  90  using  Digital  Visual  Fortran  Professional  Edition 
5.0A  and  Microsoft  Developer  Studio  97.  It  uses  NAMELIST  I/O  for  input  from  the  user 
and  produces  ASCII  text  files  that  describe  the  distribution  of  the  plasma.  From  the  input 
files  the  user  can  turn  on  or  off  certain  types  of  collisions,  customize  the  output  of  the 
program,  change  the  energy,  power,  and  other  electron  beam  characteristics,  change 
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altitude  and  velocity,  and  much  more.  Output  files  provide  descriptive  statistics  of  the 
electron  distribution  as  well  as  a  two  dimensional  map  of  the  electron  distribution. 

Design  Philosophy 

The  core  of  the  EBS  program  is  the  Monte  Carlo  method.  EBS  follows  the 
trajectory  of  one  electron  at  a  time  as  it  propagates  through  the  air.  Total  cross  sections 
from  the  theoretical  models  in  Chapter  III,  are  used  to  determine  if  the  electron 
experiences  an  elastic  or  ionization  collision.  If  an  ionization  collision  occurs  then  a  new 
electron  is  created  with  a  certain  energy  and  direction  depending  on  the  probability 
distribution  function.  As  a  result  of  the  ionization  collision,  the  incident  electron  is 
deflected  and  losses  energy  from  the  collision.  If  an  elastic  collision  occurs  then  the 
electron  is  deflected  according  to  the  probability  distribution  function  given  by  Chapter 
III,  equation  (12).  The  process  of  the  electron  colliding  with  a  neutral  molecule, 
scattering,  losing  energy,  and  moving  to  the  next  collision  is  repeated  until  the  initial 
electron’s  energy  is  less  than  the  minimum  ionization  energy  of  oxygen  (which  is  12.3 
eV)  at  which  time  the  coordinates  of  the  end  electron  are  recorded  and  a  new  electron  is 
introduced  into  the  scenario.  The  electron  termination  coordinates  and  energies  are  then 
processed  and  reported  to  the  user  via  descriptive  statistics.  For  a  flow  diagram  of  the 
EBS  program  see  Figure  23  and  24. 
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Figure  23.  Flow  Diagram  of  Main  Program. 
Numbers  represent  the  subroutines  found  in 
Table  12  that  perform  the  task  in  the  block 
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Figure  24.  Flow  Diagram  of  SimControl. 
Numbers  Represent  the  Subroutines  Found  in 
Table  12  that  Perform  the  Task  in  the  Block 
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Inputs 


EBS  is  a  data  driven  simulation  with  three  main  data  files  that  are  used  to  supply 
input  to  the  simulation.  These  files  are  formatted  for  use  with  the  Fortran  90  NAMELIST 
I/O  function  (see  Appendix  B).  The  main  data  file  is  the  Ebeam.inp  file,  which  is 
modified  by  the  user  to  provide  parameters  for  the  simulation.  The  second  data  file  is  the 
Default.dat  file,  which  provides  the  default  parameters  for  the  simulation  if  none  are 
specified  in  the  Ebeam.inp  file.  The  third  data  file  contains  data  on  the  molecules  in  the 
simulation.  Other  input  files  contain  the  ejected  electron  momentum  transfer  probability 
distribution  that  is  used  by  the  simulation  to  determine  how  much  momentum  and  energy 
is  transferred  from  an  incident  electron  to  the  ejected  or  ionized  electron. 

The  execution  of  the  Electron  Beam  Simulation  is  accomplished  by  first 
modifying  the  parameters  of  the  EBeam.inp  file  and  then  executing  the 
EbeamScattering.exe  file.  See  Appendix  B,  Table  13  for  a  listing  of  all  the  input 
variables  to  the  EBS  program  and  their  function.  Due  to  the  flexibility  of  the 
NAMELIST  I/O  format,  data  does  not  have  to  be  entered  in  a  specific  order.  If  a  number 
of  different  simulations  are  going  to  be  run,  the  user  can  modify  the  default  values  in  the 
default_file.dat  to  reduce  the  number  of  parameters  that  need  to  be  supplied  in  the 
EBeam.inp  file.  The  Molecule_Data_File.dat  can  be  modified  to  incorporate  additional 
gas  species  into  the  EBS  program.  The  data  required  to  add  a  molecule  to  the  simulation 
includes  binding  energies  of  the  electrons,  orbital  energies  of  the  electrons,  and  shell 
occupation  number.  Experimental  as  well  as  calculated  data  on  the  binding  and  orbital 
energies  of  electrons  in  various  light  molecules  and  atoms  can  be  found  on  the  NIST  web 
site  (Kim,  2000). 
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Outputs 


There  are  two  main  output  files  for  the  EBS  simulation.  The  output  files  include 
the  descriptive  statistics  file,  named  by  the  user,  which  provides  a  summary  of  all  the 
input  parameters,  the  axial  beam  profile,  the  sample  mean  of  the  transverse  electron 
distribution,  and  the  standard  deviation  of  the  transverse  electron  distribution.  The 
second  file  produced,  *. PROCESSED,  contains  a  two  dimensional  map  of  the  frequency 
distribution  of  the  coordinates  of  the  end  electrons.  Figure  25  provides  a  diagram  of  the 
coordinate  system  and  the  two  dimensional  map  in  which  the  end  electron  frequency  data 
is  stored.  The  size  and  number  of  cells  is  determined  by  the  user  and  should  be  based  on 
the  initial  energy  of  the  electron  such  that  it  provides  adequate  resolution  for  the  electron 
distribution.  Additional  files  contain  a  compilation  of  the  data  from  many  simulation 
runs  and  are  created  if  the  appropriate  options  are  selected.  This  was  done  because  many 
of  the  simulations  that  were  executed  required  too  much  computer  resources;  therefore 
large  simulations  were  divided  into  smaller  simulations.  The  results  of  the  smaller 
simulations  are  combined  into  the  *. COMPILED  and  *.RES.  The  form  of  the  data  from 
these  files  is  the  same  as  the  main  output  file  and  the  *. PROCESSED  file. 
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Axis  -  x 


•  -  Thermal  Electrons 


Figure  25.  EBS  Setup,  Electron  Beam  Propagation 
Direction,  and  Coordinate  System 


Description  of  Functions  and  Subroutines 

In  Table  10,  there  is  a  brief  synopsis  of  the  functions  and  subroutines  found  in  the 
EBS  program.  The  first  column  provides  the  number  that  is  referenced  in  Figures  23  and 
24.  The  second  column  provides  the  name  of  the  subroutine  or  function.  The  third 
column  indicates  the  main  output  of  the  function  or  subroutine  and  provides  a  brief 
description.  A  more  detailed  documentation  on  the  EBS  program  is  provided  in  the  code 
itself. 
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Table  10.  Description  of  Subroutines  for  the  EBS  Program 


Sub 

# 

Subroutine  Name 

Brief  Description 

1 

ebeamprog 

Main  Program 

The  ebeamprog  performs  four  major  functions.  First 
it  reads  the  NAMELIST  input  files  supplied  by  the 
user.  It  then  checks  the  input  data  for  errors  and 
starts  the  Monte  Carlo  simulation  by  calling  the 
simcontrol  function  if  there  are  no  errors.  When  the 
simcontrol  function  completes,  ebeamprog  outputs 
the  results  of  the  simulation  to  a  file.  Ebeamprog  also 
calls  the  functions  that  compile  the  results  of  multiple 
simulations  with  the  appropriate  input. 

EBeamSimControl 

Module 

2 

SimControl 

Monte  Carlo  Simulation 

The  controlling  function  for  the  Monte  Carlo 
simulation.  It  calls  all  the  functions  responsible  for 
setting  the  initial  energy  and  position  of  the  electrons, 
calculating  collision  cross  sections,  calculating  the 
electron  trajectory,  determining  the  number  of  ionized 
electrons  and  processing  the  results  of  the  simulation. 

EBeamFuncModule 

Module 

3 

AirDensity 

Calculate  Nm 

Calculates  the  density  of  the  air  at  the  altitude  given 
by  the  user.  AirDensity  assumes  the  atmosphere  is 
exponential  and  uses  a  scale  height  of  8180  m 
(Al’pert,  1960:84) 

4 

CreateEData 

Set  electron  initial  conditions 

Assigns  initial  energies  and  positions  to  electrons.  If 
the  logical  variable  InitDist  is  FALSE  then  the 
electrons  are  all  given  the  same  energy  and  initial 
angle.  If  the  InitDist  is  TRUE  then  the  electrons  are 
given  a  Gaussian  energy  distribution  and  initial  angle 
distribution. 

5 

GetElasticCrossSect 

Calculate  odastic 

Uses  Chapter  III,  equation  (12)  to  determine  the  total 
elastic  cross  section  of  the  electron  colliding  with 
either  nitrogen  or  oxygen  molecules. 
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GetlonCrossSect 

Calculate  cr,r n 

Uses  the  relativistic  Bethe  equation,  Chapter  III, 
equation  (20),  to  determine  the  total  ionization  cross 
section  of  an  air  molecule  based  on  the  user  supplied 
average  ionization  energy. 

7 

GetBEBCrossSect 

Calculate  alm 

Calls  the  BEBModel  function  and  passes  it  the 
appropriate  electron  orbital  data  to  calculate  the 
ionization  cross  section  of  the  orbital  shells  of  an 
atom  or  molecule.  The  orbital  shell  cross  sections  are 
then  put  into  an  array  and  passed  to  simcontrol. 

8 

SelectTime 

Calculates  time  between  collisions 

Uses  equation  (30)  to  calculate  the  time  between 
collisions  for  the  electron. 

9 

EMotion 

Calculate  e  trajectory 

Calculates  the  trajectory  of  the  electron  after  it 
undergoes  a  collision  using  the  electron  scattering 
angle  and  the  time  between  collisions  calculation. 

10 

ElasticScatteringAngle 

Calculates  elastic  scattering  angle 

Determines  the  angle  at  which  the  electron  is 
scattered  after  experiencing  an  elastic  collision.  This 
is  done  by  calling  the  function 
IntegratedElasticCrossSect,  which  uses  the  integrated 
form  of  the  Mott  elastic  cross  section  to  determine  the 
cross  section  of  an  electron  being  scattered  into  a 
particular  angular  range.  Sixty  increments  from 
thetamin  to  pi  (thetamin  is  determined  by  Chapter  III, 
equation  (13))  are  added  up  and  divided  by  the  total 
elastic  scattering  cross  section  to  determine  the 
probability  of  the  electron  being  scattered  in  that 
angular  range.  Once  the  angular  range  is  determined 
the  actual  angle  of  scattering  is  determined  by 
linearly  interpolating  between  the  two  end  values  of 
the  range  of  angles. 

11 

Ei  ectedElectronEnerg  v 

Calculate  W 

Determines  the  electron  ejection  energy  by  using  a 
previously  calculated  electron  ejection  energy 
probability  distribution  table.  The  probability 
distribution  table  was  calculated  by  integrating 

Chapter  III,  equation  (30)  over  all  ejection  angles  and 
ejection  energies.  The  results  of  the  integration  were 
stored  in  several  files  that  are  imported  into  EBS 
upon  execution  of  the  program. 
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12 

EVelocity 

Calculate  vf 

Determines  the  relativistic  electron  velocity  based  on 
the  electron  energy 

13 

IonScatteringAngle 

Calculate  e  scattering  angle 

Determines  the  angle  at  which  the  electron  is 
scattered  after  experiencing  an  ionization  collision. 

This  is  done  by  numerically  integrating  the  function 
EjectionEquation,  which  contains  Chapter  HI, 
equation  (30),  to  determine  the  cross  section  of  an 
electron  being  scattered  into  a  particular  angular 
range.  Sixty  increments  from  0  to  thetamax 
(thetamax  is  set  such  that  the  probability  of  the 
electron  being  scattered  at  an  angle  larger  than 
thetamax  is  10'7)  are  added  up  and  divided  by  the 
total  scattering  cross  section  to  determine  the 
probability  of  the  electron  being  scattered  in  a 
particular  angular  range.  Once  the  angular  range  is 
determined  the  actual  angle  of  scattering  is 
determined  by  linearly  interpolating  between  the  two 
end  values  of  the  range  of  angles. 

14 

Ej  ectedElectron  Angle 

Calculate  e  scattering  angle 

Same  as  IonScatteringAngle,  except  the 

MottEquation  function  is  used  which  contains 

Chapter  HI,  equation  (23). 

15 

AngleCleanUp 

Keens  angles  between  180  and -180  degrees 

Alters  the  scattering  angle  of  the  electron  if  it  is 
greater  than  180  degrees  or  less  than  -180  degrees 
such  that  it  is  in  the  appropriate  quadrant  and  has  a 
value  of  180  to  -180  degrees. 

16 

! 

MaxE 

Determines  maximum  number  e  simulated 

Determines  the  maximum  number  of  electrons  that 
will  be  modeled  by  the  Monte  Carlo  simulation  by 
dividing  initial  electron  energy  by  the  minimum 
ionization  energy 

FunctionModule 

Module 

17 

ScatteringAngleDist 

Equation 

Contains  Chapter  HI,  equation  (12)  which  describes 
the  elastic  scattering  angle  distribution  of  an  electron 
after  colliding  with  a  nucleus 

18 

IntegratedElastic 

CrossSect 

Equation 

Contains  Chapter  III,  equation  (14)  which  describes 
the  elastic  scattering  angle  distribution  of  an  electron 
after  colliding  with  a  nucleus 
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19 

EjectionEquation 

Equation 

Contains  Chapter  III,  equation  (23)  which  is  Mott’s 
ionization  TDCS 

20 

MottEquation 

Equation 

contains  Chapter  III,  equation  (36)  which  is  Mott’s 
ionization  DDCS 

21 

BEBModel 

Equation 

Contains  Chapter  III,  equation  (48)  which  is  the 
RBEB  total  ionization  cross  section 

IntegrationModule 

Module 

22 

Integrate 

Integrates  equations 

Integrates  functions  using,  Simpson’s,  Trapezoid, 
or  Gaussian  Quadrature  method. 

EBeamlnput 

Module 

23 

GetCrossSectData 

Imports  cross  section  data 

Imports  normalized  ejection  energy  probability 
distribution  data 

24 

CheckDatal 

Checks  data  for  errors 

Checks  data  for  errors  and  sets  error  flag  if  error  in 
input  data  is  found 

25 

CheckData2 

Create  error  message 

Generates  appropriate  error  message  if  error  in 
input  data  is  found 

EBeamOutput 

Module 

26 

ProcessOutput2 

Process  Monte  Carlo  results 

Takes  the  analog  electron  distribution  calculated  by 
the  Monte  Carlo  simulation  and  places  it  in  a  two 
dimensional  grid  of  the  electron  beam  based  on  the 
end  x  and  y  coordinates  of  the  electrons. 
ProcessOutput2  then  determines  the  axial  profile, 
the  mean  transverse  electron  position,  and  the 
transverse  electron  distribution. 

27 

ProcessCompiledOutput 

Process  Monte  Carlo  results 

Same  as  ProcessOutput2  except  used  to  compile 
results  from  multiple  simulations 

28 

ReCompileOutput 

Process  Results 

Combines  archived  data 

A  validation  of  the  Electron  Beam  Scattering  (EBS)  simulation  was  performed  to  insure 
that  the  Monte  Carlo  method  was  implemented  correctly  in  the  EBS  program.  The  EBS 
simulation  was  compared  with  the  results  of  the  SEBPM.  For  the  purposes  of  this 
comparison  the  electrons  in  the  EBS  simulation  were  limited  to  the  same  restrictions 
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applied  to  the  RNE  scenario  of  the  SEBPM  model;  that  is  the  electrons  did  not 
experience  any  angular  scattering  and  the  ejected  electrons  received  no  energy.  In  Figure 
26,  a  comparison  between  the  electron  distribution  profile  of  the  SEBPM  and  EBS 
programs  is  shown.  The  results  of  the  simulations  showed  excellent  correlation 
indicating  that  the  cross  section,  time,  and  trajectory  calculations  in  the  Monte  Carlo 
simulation  had  been  implemented  correctly.  The  results  also  indicated  that  the  SEBPM 
calculations  were  indeed  the  bounding  conditions  for  the  electron  beam  longitudinal 
extent  and  electron  number  profile. 


Figure  26.  Comparison  between  the  SEBPM  and  EBS  Simulations 
Axial  Electron  Density  Profile  for  the  RNE  scenario 


The  EBS  program  can  provide  a  relatively  accurate  estimate  of  the  plasma  density 
distribution  and  the  dimensions  of  the  electron  beam  generated  plasma.  By  importing  the 
results  of  the  EBS  program  in  to  the  EMWPM  program,  we  can  obtain  power  attenuation 
due  to  a  plasma  that  has  achieved  a  steady  state  in  density.  However,  the  EBS  simulation 
completely  excludes  any  volumetric  loss  mechanisms  introduced  by  reactions  between 
the  electrons,  ions,  and  neutral  molecules  present  in  the  plasma.  The  next  section 
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provides  an  estimate  as  to  the  plasma  density  variations  due  to  attachment,  detachment, 
and  recombination  processes. 

Plasma  Density  Loss  Mechanisms 

Now  that  we  have  the  means  to  calculate  the  source  term  for  the  electrons  being 
injected  into  the  plasma,  we  must  consider  the  plasma  chemistry  that  results  from  having 
a  highly  reactive  ion  and  electron  mix.  The  following  section  describes  the  major  loss 
and  gain  mechanisms  in  a  nitrogen-oxygen  plasma  and  the  resulting  temporal  evolution 
in  the  concentration  of  the  electrons,  ions,  and  neutral  atoms.  The  chemical  rate 
coefficients  associated  with  the  relevant  kinetic  mechanisms  are  assembled  into  a  set  of 
differential  equations  that  describe  the  rate  of  change  in  the  density  of  the  constituents  of 
the  plasma.  From  these  differential  equations,  we  obtain  an  estimate  of  the  variation  in 
the  density  of  the  plasma  as  a  function  of  time.  For  simplification  a  large  number  of 
insignificant  reaction  processes  were  excluded  from  the  model.  Also  atoms  or  molecules 
in  an  excited  state  after  a  reaction  were  assumed  to  de-excite  and  join  the  ground  state 
population  of  atoms  or  molecules  immediately. 

As  the  electrons  decrease  in  energy,  collisions  other  than  ionization  and  elastic 
scattering  become  relevant  to  our  calculation  of  the  plasma  density.  For  electrons  the 
loss  and  gain  mechanisms  other  than  ionization  are  attachment,  recombination,  and 
detachment.  Attachment  occurs  when  a  free  electron  becomes  bound  to  a  neutral  atom 
or  molecule,  forming  a  negative  ion.  The  following  equation  describes  an  electron 
attaching  to  an  oxygen  molecule  in  a  two-body  attachment  process. 

6  +  02  — ^  C?2  (47) 
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At  higher  pressures  of  02 ,  another  attachment  process,  a  three-body  attachment  process, 
may  dominate: 

e  +  2 02->0~+02  (48) 

Recombination  describes  the  process  of  an  electron  colliding  with  a  positive  ion  and 
forming  a  neutral  atom  or  molecule 

6  +  02  — ^  02  (49) 

Detachment  describes  the  process  of  an  electron,  atom  or  molecule  colliding  with  a 
negative  ion  and  stripping  the  attached  electron. 

<?  +  02  — ^  02  +  2e  (50) 

The  dominant  loss  mechanisms  for  the  plasma  based  on  their  large  rate  constants 
and  the  density  of  the  reactants  are  the  attachment  processes.  The  following  attachment 
processes  are  the  most  dominant  of  all  the  attachment  processes 


c  +  02  — ^  O  0 

(51) 

€  +  202  — ^  @2  ^2 

(52) 

The  first  reaction  is  a  dissociative  attachment  and  the  second  reaction  is  three-body 
molecular  attachment.  Dissociative  attachment  is  most  prevalent  at  electron  energies 
between  4  to  12  eV,  which  is  right  below  the  primary  ionization  energy  of  molecular 
oxygen.  The  three-body  molecular  attachment.  Equation  (48),  is  most  prevalent  at 
electron  energies  of  0.1  to  1  eV.  Molecular  and  atomic  nitrogen  do  not  form  stable 
negative  ions,  therefore  nitrogen  attachment  rates  are  negligible  and  not  considered  in  the 
calculations.  Recombination  reactions  occur  predominantly  at  electron  energies  of  less 
than  0.1  eV  or  temperatures  of  less  than  910  K  (See  Appendix  C).  Some  molecular 
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detachment  reactions  are  very  temperature  dependent  (such  as  reactions  14  and  15  in 
Appendix  C),  hence  the  rate  coefficient  of  the  reaction  increases  dramatically  with 
increased  temperature.  The  electron  affinity  between  the  attached  electron  and  atomic 
and  molecular  oxygen  is  1.465  eV  and  0.44  eV  respectively.  As  a  result  the  electron 
detachment  processes  only  occurs  if  the  incident  electron  energy  is  greater  than  the 
attached  electron  affinity  to  atomic  or  molecular  oxygen.  Since  the  concentration  of  both 
positive  and  negative  ions  is  very  important  for  determining  the  recombination  and 
detachment  rates  in  the  plasma,  the  ion-ion  and  ion-neutral  reactions  must  be  considered 
as  well.  Most  of  the  reactions  in  the  calculations  are  not  described  in  this  section  for  the 
sake  of  brevity,  but  a  complete  listing  of  all  the  reactions  used  in  the  plasma  chemistry 
calculations  can  be  found  in  Appendix  C. 

The  reactions  in  Appendix  C  were  selected  to  be  in  the  model  because  they  would 
have  a  significant  impact  on  the  density  of  the  constituents  of  the  plasma.  In  general,  any 
reaction  between  molecular  oxygen  or  nitrogen  and  an  ion  or  electron  was  selected  to  be 
a  reaction  in  the  model  because  of  the  high  density  of  molecular  oxygen  and  nitrogen. 
Atom-ion  and  ion-ion  reactions  were  included  if  the  rate  constant  of  the  reaction  was 
sufficiently  large  that  the  reaction  may  have  a  noticeable  effect  on  the  concentrations  in 
the  plasma. 

From  the  reactions  in  Appendix  C,  17  single  order  differential  equations  were 
generated  to  calculate  the  change  in  the  concentrations  of  the  ions  and  electrons  in  the 
plasma.  To  better  model  the  electron  attachment  and  detachment  processes  that  occur  in 
certain  energy  ranges,  the  electrons  were  divided  into  three  energy  groups  of  high, 
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medium,  and  low  energies.  Where  the  low  energy  is  below  0.1  eY,  medium  is  between 
0.1  and  4  eV,  and  high  is  between  4  eV  and  12.3  eV  (see  Table  11  for  more  details). 


Table  11.  Electron  Energy  Ranges 


Classification 

Energy  Range 

Reactions 

High 

4  to  12.3  eV 

e**  +  02  —>0  +  0~ 
e**  +  02  —>  02  +2e 

e**  +  0~  ~^0  +  2e 

Medium 

0.1  to  4  eV 

6  *  +2 02  — ^  02  +  02 

6  *  +(?2  +  N 2  — ^  N2  +  02 
e*+02~  ->02  +  2e 

Low 

<0.1  eV 

See  Appendix  C 

The  rate  equations  were  developed  from  the  reactions  in  Table  1 1  by  treating  the 
electrons  in  different  energy  ranges  as  different  elements.  This  was  done  so  that  the 
population  of  electrons  in  a  particular  energy  range  could  only  react  with  the  dominant 
processes  of  that  energy  range.  The  high  and  medium  energy  ranges  consist  of  the 
dominant  attachment  and  detachment  process  for  that  energy  range.  If  a  high  energy 
electron  detaches  an  electron  from  an  ion,  then  the  ejected  and  scattered  electrons  have  a 
probability  of  being  in  any  of  the  three  energy  ranges  based  on  Mott’s  ejection  energy 
probability  distribution  (Chapter  III,  equation  (30)).  If  a  medium  energy  electron 
detaches  an  electron  from  an  ion,  then  the  ejected  and  scattered  electrons  have  a 
probability  of  being  in  either  the  medium  or  the  low  energy  range  based  on  Mott’s 
ejection  energy  distribution. 

The  method  of  dividing  the  electron  energies  into  different  categories  was 
required  to  obtain  a  better  estimate  of  the  electron  densities  in  the  plasma.  This  is  because 
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reaction  1  in  Appendix  C  is  only  significant  if  the  electrons  have  a  high  energy  (greater 
than  4  eV).  Since,  only  a  few  percent  of  the  electrons  have  sufficient  energy  to  react  in 
reaction  1,  electrons  at  higher  energies  were  separated  from  electrons  at  lower  energies  so 
that  only  they  would  be  involved  in  reaction  1.  The  rate  constants  for  reactions  3  and  4 
are  very  different  for  electrons  at  low  energies  versus  electrons  at  higher  energies. 
Therefore,  reactions  3  and  4  were  included  in  both  the  medium  and  low  energy  ranges 
with  a  rate  constant  that  was  appropriate  for  the  electron  energy  range. 

If  we  assume  that  the  production  rate  of  thermal  electrons  into  any  point  in  the 
plasma  density  is  constant,  then  the  electron  beam  can  be  treated  as  a  constant  source  of 
new  electrons  into  the  reaction.  The  electron  beam  also  provides  a  source  term  for  the 
positive  oxygen  and  nitrogen  ion  rate  equations  because  an  equal  number  of  positive  ions 
and  electrons  must  be  made  (assuming  that  single  ionization  dominates).  However,  it 
was  unknown  how  many  of  the  thermal  electrons  will  fall  within  the  energy  ranges  in 
Table  11.  To  approximate  the  thermal  electron  distribution,  the  Monte  Carlo  simulation 
was  run  to  obtain  an  approximate  end  electron  energy  distribution.  The  resulting  electron 
energy  distribution  was  exponential  and  when  integrated  indicated  that  approximately  1/3 
of  the  electrons  fell  within  each  energy  range  listed  in  Table  11.  Therefore,  the  constant 
source  term  for  each  of  the  electron  rate  equations  is  given  by  the  following  expression 


(53) 


where 


y123  =  electron  density  source  term  for  a  particular  energy  range 


y0  =  rate  of  electron  density  flow  into  a  volume 
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At  =  pulse  time  length 

No  cross  section  data  was  available  for  the  electron  detachment  processes 
described  in  Table  11,  therefore  a  modified  BEB  model  was  used  to  estimate  the 
detachment  cross  section  for  0{  and  O'.  According  to  Kim,  better  agreement  between 
the  BEB  model  and  experimental  data  on  the  ionization  of  ions  was  obtained  when  the  ad 
hoc  term  in  the  denominator  of  Chapter  III,  equation  (42)  was  changed  from  (T  +  B  +  U) 
to  T  +  (B  +  U)/2  at  non-relativistic  T  (Kim,  2000:052710-5).  The  BEB  model  yielded  the 
detachment  cross  section  of  O'  and  Oi  depicted  in  Figure  27. a  and  27.b  respectively.  By 
obtaining  the  electron  detachment  cross  section,  a  rate  coefficient  for  the  electron 
detachment  process  can  be  calculated  using  the  relationship 

k  =  ov  (54) 


where 


k  =  reaction  rate  constant 


v  =  the  velocity  of  the  electrons 


the  average  rate  constant  over  the  energy  ranges  of  interest  was  determined  for  the 
reactions  in  Table  11. 
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Figure  27.  Detachment  Cross  Sections 
for  a)  O'  b)  Oi 
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The  rate  equations  developed  from  Table  1 1  were  combined  with  the  rate 
equations  developed  from  the  chemical  reactions  in  Appendix  C  to  model  the  temporal 
evolution  of  the  plasma.  Due  to  the  enormous  number  of  reactions  that  are  modeled,  a 
Mathematica  code  was  developed  that  converted  the  chemical  reactions  in  Appendix  C 
into  the  17  single  order  differential  rate  equations  that  are  used  to  model  the  plasma 
chemistry.  The  resulting  rate  equations  were  solved  using  NDSolve,  the  numerical 
differential  equation  solver  in  the  software  package  Mathematica. 

Results  of  Plasma  Chemical  Reaction  Calculations 

The  results  of  the  plasma  chemistry  model,  which  was  developed  in  the  previous 
section  for  various  electron  and  ion  density  source  terms,  are  shown  in  Appendix  D.  In 
Figure  39,  the  resulting  electron  and  ion  densities  due  to  the  electron  beam  being  on  for  5 
ms  and  then  shutting  off  are  shown.  The  results  of  these  calculations  were  used  to 
modify  the  plasma  density  obtained  from  the  EBS  program.  This  is  done  so  that  a  better 
estimate  can  be  obtained  for  the  attenuation  and  refraction  due  to  a  plasma  made  by  an 
electron  beam  with  a  certain  power  and  initial  electron  energies. 

From  the  figures  in  Appendix  D,  we  see  that  the  electron  density  reaches  an 
approximate  steady  state  in  a  few  microseconds  regardless  of  the  magnitude  of  the 
electron  density  source  terms.  The  steady  state  is  due  mostly  to  a  balance  between  the 
three-body  attachment  reactions  3  and  4  and  the  electrons  generated  by  the  electron 
beam.  From  Figures  40.a,  40.b,  and  40.c,  we  see  that  the  steady  state  density  of  the 
electrons  is  directly  proportional  to  the  magnitude  of  the  source  term  of  electrons  from 
the  electron  beam  into  the  reaction.  This  relationship  can  easily  be  seen  if  we  consider 
that  the  loss  of  electron  density  due  to  attachment  is  much  greater  than  any  other  loss 
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mechanisms  in  the  plasma  and  that  the  electron  density  reaches  a  steady  state.  From  these 
assumptions,  we  can  approximate  the  time  rate  of  change  of  the  electron  density  as 

=  -«JV202  -  kalN,No2NK2  +  y,  =  0  (55) 

at 

Ne  = - — - =  — — —  (56) 

N  02^3  +  Nq2N N2^4  ^ attached 

where 

Ne  =  electron  density 
N02  =  molecular  oxygen  density 
NN2=  molecular  nitrogen  density 
yl  =  electron  source  term  due  to  the  electron  beam 
ki  =  rate  constant  for  reaction  3 
kA  =  rate  constant  for  reaction  4 
^attached  =  rate  electrons  attach  to  molecular  oxygen 
The  terms  in  the  denominator  of  equation  (56)  are  all  constant;  therefore  the  electron 
density  is  directly  proportional  to  the  electron  beam  source  term,  which  is  consistent  with 
the  data  presented  in  Figure  39  and  Figure  28.  This  relationship  holds  for  y1  values  from 

10n  to  10 22  (e~  / s -cm3). 
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Figure  28.  Log  of  Steady  State  Electron  Density 
Versus  Log  of  Electron  Beam  Source  Term  y 

Calculations  for  electron  densities  at  various  pulse  durations  were  used  to 
determine  an  exact  relation  between  the  rate  at  which  the  electron  beam  injects  electrons 
and  ions  into  the  plasma  and  the  steady  state  electron  density.  In  Figure  29,  the  y-axis  is 
the  log  of  the  ratio  of  the  steady  state  electron  density,  Ne,  with  plasma  chemistry  losses 
to  the  plasma  density  achieved  by  pulsing  the  electron  beam  for  a  period  of  tp ,  resulting 

in  an  electron  density  of  yrp  for  the  case  of  no  loss  mechanisms.  The  reason  for  this 

analysis  is  to  provide  a  simple  way  of  translating  the  electron  beam  density  distribution 
without  losses  into  a  density  distribution,  which  includes  losses  due  to  plasma  chemistry. 
This  treatment  is  only  appropriate  because  the  electron  density  with  losses  achieves  a 
near  steady  state  in  a  few  microseconds;  therefore  the  electron  density  with  losses  is 
approximately  constant  for  the  duration  of  the  pulse. 

The  results  of  those  calculations  over  the  range  of  pulse  lengths  between  0.0001  s 
and  0.005  s  are  shown  in  Figure  29.  We  see  that  the  shorter  the  pulse  duration  (therefore. 


108 


a  higher  rate  of  electrons  being  inserted  into  the  plasma),  the  lower  the  steady  state 
electron  density  losses  due  to  attachment  and  recombination.  From  Figure  29,  we  see 
that  the  percent  loss  of  electron  density  is  approximately  the  same  for  all  electron  and  ion 
injection  rates  from  the  electron  beam,  therefore,  we  can  approximate  the  plasma  loss 
mechanisms  as  having  a  constant  percent  loss  for  all  electron  densities  in  the  plasma.  As 
a  result,  a  constant  Plasma  Density  Loss  Factor  (PDLF)  can  be  applied  to  the  plasma 
density  without  losses  to  approximate  a  plasma  density  with  losses. 

It  should  be  noted  that  there  is  a  very  large  overall  reduction  in  the  electron 
density  of  the  plasma  due  to  the  attachment  and  recombination  reactions.  For  the 
particular  electron  and  ion  density  source  terms  shown  in  Figure  29.a,  the  ratio  of  the 
steady  state  plasma  density  with  and  without  losses  is 

—  =  -^-  =  3xl0~5  (57) 

N0  YSp 

where 

N0  =  the  plasma  density  without  losses  after  a  time  rp 
Tp  =  pulse  width 

If  we  assume  again  that  all  other  loss  or  gain  mechanisms  are  negligible  compared  to  the 
three  body  attachment  processes  then,  then  we  may  substitute  equation  (56)  into  equation 
(57)  to  obtain  the  expression 

^l  = - £ - = - £ -  (58) 

Nq  (N  02&3  +  N q2N #2^4 ^  attached  ^  p 

Since,  Rattached  is  approximately  constant  then  the  ratio  of  the  electron  densities  with  and 
without  losses  scales  as  the  inverse  of  the  pulse  time.  From  Figure  29,  we  can  tell  that 
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this  relationship  holds  exactly  for  lower  electron  injection.  The  large  loss  in  electron 
densities  is  due  to  molecular  oxygen  and  nitrogen  being  in  very  plentiful  supply.  As  a 
result  the  majority  of  electrons  become  attached  to  molecular  oxygen,  hence  reducing  the 
electron  density  by  several  orders  of  magnitude  through  out  the  plasma.  Also  as  a  result 
of  the  rapid  attachment  rate,  the  plasma  becomes  a  positive  and  negative  oxygen  ion 
plasma,  which  may  result  in  the  negative  and  positive  oxygen  ions  becoming  the 
dominant  term  in  determining  the  plasma  frequency. 

In  Figure  29,  the  slight  decrease  in  the  ratio  of  electron  densities,  Ne/N0  ,  at 
medium  injection  rates  is  due  to  recombination  reactions  7, 10,  and  46,  which  contribute 
to  the  electron  losses  due  to  the  higher  densities  of  positive  ions.  The  ratio  Ne/N0 
increases  rapidly  at  the  highest  electron  injection  rates  because  the  electron  densities 
increase  fast  enough  that  the  detachment  processes  can  free  a  significant  portion  of  the 
electrons  bound  to  the  negative  oxygen  ions. 


a)  b) 
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Figure  29.  Ratio  of  Electron  Densities  with  and  without  Loss  Mechanisms 
versus  the  Log  of  the  Initial  Electron  Density  at  Different  Pulse  Lengths  of 
a)  0.005  s  b)  0.0025  s  c)  0.0005  s  d)  0.0001  s 


It  should  also  be  noted  from  the  figures  in  Appendix  D  that  the  electron  densities 
diminish  very  rapidly  when  the  electron  beam  is  turned  off.  However,  from  Figure  39, 
we  see  that  the  negative  molecular  oxygen  ion  density  remains  fairly  constant  after  the 
electron  beam  is  shut  off.  Since  the  electron  densities  are  negligible,  the  plasma 
frequency  will  be  a  function  of  the  negative  oxygen  ion  density  whose  plasma  frequency 
is  approximately  two  orders  of  magnitude  less  than  the  plasma  frequency  for  an 
equivalent  electron-positive  ion  plasma.  To  achieve  a  negative  molecular  oxygen  ion 
plasma  frequency  near  a  GHz,  and  hence  a  significant  degree  of  attenuation  of  an  EM 
wave  in  the  GHz  range,  after  the  electron  beam  has  been  turned  off  requires  an 
extraordinarily  high  negative  ion  density  (1013-1014  e/cm3).  Therefore,  when  the  electron 
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beam  is  off,  the  plasma  density  will  not  be  sufficient  to  attenuate  or  refract  an  EM  wave 
in  the  range  of  frequencies  presented  in  this  study. 

The  calculations  of  the  last  section  were  a  rough  estimate  of  the  reactions  that 
would  occur  in  the  plasma.  The  air  molecules  are  assumed  to  be  unheated  and  have  a 
temperature  of  300  K,  however  this  assumption  may  not  hold  if  the  electron  beam  is 
operated  for  a  long  period  of  time  at  a  higher  power  setting.  The  change  in  electron 
energies  was  also  modeled  very  coarsely  with  the  high,  medium,  and  low  energy  ranges 
and  should  be  examined  rigorously  with  a  Boltzmann  transport  calculation  to  obtain  a 
better  estimate  of  the  temporal  evolution  of  the  plasma  density. 

In  this  last  section,  we  have  developed  a  means  to  gauge  the  loss  of  plasma 
density  due  to  attachment  and  recombination  mechanisms.  The  ratio  of  the  electron 
densities  of  the  plasma  with  and  without  losses  was  shown  to  be  constant  over  the  entire 
density  range  of  the  plasma,  therefore  an  attenuation  factor  can  be  applied  uniformly  to 
the  plasma  density  without  losses  to  estimate  the  effects  of  the  loss  mechanisms  on  the 
plasma  density. 
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V.  Results  and  Conclusions 


This  chapter  presents  results  derived  from  the  computer  programs  Electro- 
Magnetic  Wave  Propagation  Model  (EMWPM)  and  the  Electron  Beam  Scattering  (EBS) 
simulation  and  the  plasma  chemistry  model  described  in  the  previous  section.  The 
primary  purpose  of  this  chapter  is  to  demonstrate  the  capability  of  the  computer  models 
developed  in  the  previous  chapters  to  quantify  the  plasma  density  distribution  resulting 
from  the  injection  of  a  relativistic  electron  beam  into  air.  The  data  presented  in  this 
chapter  is  limited,  because  the  parameter  space  for  the  scenario  is  quite  extensive  and 
very  dependent  on  the  environment  in  which  the  electron  beam  will  be  utilized. 

Plasma  Density  and  Spatial  Distribution 

This  section  presents  the  results  obtained  from  the  EBS  computer  program 
described  in  Chapter  IV.  The  results  of  the  density  distribution  of  the  plasma  without 
losses  mechanisms  as  predicted  by  the  EBS  simulation  are  shown.  The  results  of  the 
EBS  simulation,  which  describes  the  spatial  extent  of  the  plasma  when  using  the  electron 
beam  at  various  initial  electron  energies  and  altitudes  is  also  reported.  Finally,  the 
section  concludes  with  a  summary  of  how  the  descriptive  statistics  of  the  plasma  reported 
by  EBS  are  used  to  generate  a  plasma  density  distribution. 

Static  Plasma  Distribution 

In  Table  12,  the  run  matrix  for  the  EBS  program  at  different  altitudes  and  at 
different  initial  electron  energies  is  shown.  This  run  matrix  was  executed  in  order  to 
explore  the  parameter  space  of  longitudinal  extent,  transverse  extent,  and  densities  of  the 
electron  beam  generated  plasma.  The  main  output  of  EBS  includes  the  number  of 
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electrons  per  cm  in  the  axial  direction,  mean  and  standard  deviation  of  the  transverse 
distribution,  and  a  file  containing  a  two  dimensional  map  of  the  plasma  electron  densities. 
A  sample  of  the  results  from  two  different  EBS  simulations  at  500  keV  and  300  keV 
initial  electron  energies  are  shown  in  Figure  30.  The  simulation  parameters  for  all  the 


runs  performed  in  this  section  can  be  found  in  Table  13. 

Table  12.  Run  Matrix 


Electron 

Energy 

Altitudes 

Air  Number 
Density 

100  keV 

0m 

2.69xl019  cm'3 

300  keV 

2500  m 

1.98xl019cm'3 

500  keV 

5000  m 

1.46xl019cm'3 

750  keV 

7500  m 

1.07xl019cm'3 

1  MeV 

10000  m 

7.91xl018cm'3 

Table  13.  EBS  Parameters  for  all  Simulation  Runs 


Parameter 

Value 

Parameter 

Value 

Power 

10,000  W 

Use  Exponential 
Atmosphere 

True 

Pulse  Length 

0.5  ms 

Elastic  Collisions 

True 

Num  Simulated  Electrons 
Move  Created  Electrons 

1500 

False 

Inelastic  Collisions 

True 
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a) 

$  of  Electrons 
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#  of  Electrons 


Figure  30.  Plasma  Distribution  Results  from  the  EBS  simulation 
a)  Axial  Density  Profile  for  500  keV  e'  b)  Mean  of  Transverse 
Distribution  for  500  keY  e'  c)  STD  of  Transverse  Distribution  for 
500  keV  e' d)  Axial  Density  Profile  for  300  keV  e‘  e)  Mean  of 
Transverse  Distribution  for  300  keV  e'  f)  Mean  of  Transverse  for 
300  keV  e' 

In  Figure  30,  the  axial  and  transverse  distribution  plots  are  both  fit  with  a  fourth 
order  polynomial  and  the  mean  transverse  distribution  plot  was  fit  with  a  linear  function. 
(Which  was  the  case  for  all  the  results  from  the  EBS  simulation).  The  axial  density 
profile  predicted  by  the  EBS  simulation  is  shown  in  Figure  30.a  and  Figure  30.d.  The 
results  of  the  EBS  simulation  with  lateral  scattering  of  electrons  indicate  that  the 
maximum  in  the  axial  density  will  be  located  near  the  center  of  the  profile,  where  as  the 
results  without  lateral  scattering  show  that  the  maximum  electron  density  is  located  near 
the  end  of  the  profile,  which  are  shown  in  Figure  26.  The  peak  value  in  the  non-lateral 
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scattering  case  was  due  to  the  ionization  cross  section  achieving  a  maximum  value  at  low 
electron  energy  levels.  For  the  lateral  scattering  case,  the  few  electrons  that  travel  the 
farthest  in  the  x  direction  are  those  electrons  that  have  experienced  the  least  deflection 
and  hence  least  energy  loss  due  to  inelastic  collisions.  However,  the  majority  of  the 
electrons  experience  a  larger  degree  of  deflection  at  an  earlier  time  in  their  propagation. 
This  results  in  the  electrons  traveling  in  the  transverse  direction  to  the  beam.  Even 
though  the  distance  traveled  through  the  air  is  the  same  for  the  scattered  electrons,  they 
do  not  travel  as  far  in  the  axial  direction.  This  behavior  results  in  a  maximum  in  the 
center  rather  than  at  the  end  of  the  axial  electron  density  profile. 

In  Figure  30.c  and  30.f,  a  standard  deviation  (STD)  is  used  to  measure  the 
transverse  extent  of  the  plasma  because  the  transverse  distribution  of  the  plasma  is 
Gaussian.  The  STD  describes  the  radial  position  at  which  the  Gaussian  distribution  is  at 
half  its  maximum  value,  hence  it  is  related  to  the  transverse  extent  of  the  plasma.  The 
reason  that  the  transverse  distribution  is  Gaussian  is  described  later  on  in  this  section. 
From  Figure  30.a  and  30.d,  we  notice  that  the  length  of  the  500  keV  profile  is  nearly 
twice  the  length  of  300  keV  profile.  As  a  result,  the  number  of  electrons  per  cm  doubles 
from  Figures  30.a  to  30.d  for  the  fixed  power  condition.  We  also  notice  that  the 
transverse  standard  deviation  of  the  500  keV  electron  beam  is  roughly  twice  the 
transverse  standard  deviation  of  the  300  keV  electron  beam.  The  reason  for  the  doubling 
of  the  transverse  and  longitudinal  extent  of  plasma  is  because  the  electrons  with  500  keV 
energy  can  experience  nearly  twice  as  many  collisions  as  the  300  keV  electrons  before 
becoming  thermal  electrons,  hence  allowing  the  500  keV  electrons  to  travel  nearly  twice 
as  far.  Also  the  higher  energy  electrons  have  a  smaller  collision  cross  section,  hence  they 
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have  a  longer  mean  free  path,  resulting  in  an  overall  larger  longitudinal  and  transverse 
extent  to  the  plasma  density  distribution.  Therefore,  at  the  same  power  settings  and  pulse 
duration  the  density  of  the  300  keV  plasma  will  be  roughly  8  times  greater  than  the  500 
keV  plasma,  potentially  resulting  in  a  huge  difference  in  the  attenuation  and  refraction 
resulting  from  the  two  plasmas.  It  should  also  be  noted  that  the  shape  of  the  longitudinal 
and  transverse  electron  density  profile  is  only  a  function  of  the  initial  energy  of  the 
electrons  in  the  electron  beam.  Therefore,  if  the  initial  electron  energy  is  kept  the  same, 
then  the  electron  density  is  linearly  proportional  to  the  electron  beam  power. 

The  results  of  the  run  matrix  in  Table  12  are  shown  in  Figure  3 1 .  The  fit  to  the 
data  in  Figure  31  is  a  quadratic  function,  because  the  ionization  cross  section  decreases 
slightly  with  increased  electron  energy.  This  slight  decrease  in  cross  section  results  in  a 
slightly  longer  mean  free  path  for  higher  energy  electrons,  hence  the  higher  energy 
beam’s  longitudinal  extent  is  not  linearly  proportional  to  the  energy  of  the  electrons. 


- 5000  m 

- Sea  Level 
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Figure  31.  Maximum  Length  and  Displacement  of 
Plasma  versus  Electron  Energy  at  Various  Altitudes 
a)  Plasma  Length  b)  Plasma  Transverse  STD 


It  should  be  noted  that  Figure  3  Lb  represents  one  standard  deviation  (STD)  in  the 
transverse  plasma  distribution  (i.e.  68%  of  the  electrons  are  within  that  radius,  98%  of  the 
electrons  are  within  twice  the  radius),  therefore  there  is  still  a  non-zero  plasma  density  at 
transverse  distances  exceeding  the  values  given  in  Figure  3 l.b.  If  the  plasma  density  and 
hence  the  plasma  frequency  is  high,  then  the  plasma  density  out  to  the  second  standard 
deviation  may  be  adequate  to  attenuate  an  EM  wave. 

Figure  32  shows  how  the  standard  deviation  and  length  of  the  plasma  vary  with 
altitude.  The  function  used  to  fit  the  data  presented  in  Figure  32  is  the  inverse  of  the 
function  used  to  determine  the  number  density  of  the  air  at  various  altitudes  (exp(y  /  H) ), 
where  y  is  the  altitude  and  H  is  the  scale  height  of  the  atmosphere.  From  this  result,  we 
conclude  that  the  longitudinal  extent  and  STD  in  the  transverse  direction  are  inversely 
proportional  to  the  number  density  of  the  atmosphere.  This  is  not  surprising  because  the 
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function  controlling  the  time  that  the  electron  travels  between  collisions  (Chapter  IV, 
equation  (30))  is  inversely  proportional  to  the  number  density  of  the  plasma. 


-  300  keV 

Figure  32.  Plasma  Profile  versus  Altitude 
a)  Plasma  Length  b)  Plasma  Transverse  STD 

The  functions  obtained  by  fitting  the  data  in  Figure  30  were  used  to  develop  a 
density  map  of  the  plasma.  The  reason  for  using  the  fit  of  the  data  rather  than  the  data 
itself  is  that  the  spikes  in  the  EBS  data  create  problems  in  the  propagation  model.  Fermi 
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and  Bethe  both  predicted  that  the  transverse  distribution  of  a  beam  of  electrons 
undergoing  multiple  collisions  would  be  Gaussian  (Orear,  1950:36).  In  Figure  33, 
samples  of  the  transverse  distribution  predicted  by  EBS  are  shown. 


a)  b) 


Figure  33.  Transverse  Distribution  of  the  Electron  Beam  Generated  Plasma 
for  Electron  Energies  of  500  keV  at  Axial  Distances  of  a)  35  cm  b)  55  cm 


From  this  information  and  the  functions  fit  to  the  data  in  Figure  30,  the  following 
empirical  model  describing  the  number  of  electrons  in  the  plasma  at  a  point  was 
developed 


N(x,y)  = 


m 

^27lg(x)2 


exp[ 


2ng(x) 


(1) 


where 

N(x,  y )  =  number  of  electrons  at  coordinate  (x,y) 

/( x)  =  fit  of  the  axial  electron  number  data 
g(x)  =  fit  of  the  transverse  standard  deviation  data 
The  EBS  simulation  was  limited  to  two-dimensions,  therefore  to  calculate  the  density  of 
the  electron  beam  we  will  now  consider  the  third  dimension.  To  do  this  we  will  use  a 
cylindrical  coordinate  system  in  which  the  Cartesian  coordinates  (x,  y,z)  are  transformed 
to  cylindrical  coordinates  ( z,r,6 ) .  Since,  Mott’s  elastic  scattering  and  ionization  cross 
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sections  were  isotropic  in  the  scattering  and  ejection  angle  distribution  of  the  angles  <j> 


and  y/  respectively  (see  Chapter  III,  Mott’s  Ionization  Cross  Section),  the  electron  beam 
will  be  isotropic  in  the  angle#  as  well.  To  obtain  the  plasma  density  for  a  cell  in  the 
plasma,  the  following  equation  was  used 


P  = 


Z2r2 

JJ 

zl  rl 


fU)  -exp[-  ^ 


V2fl-g(z)2  2;rg(z): 

z2r 2 

—  j  J  Imdrdz 


-]drdz 


(2) 


zl  rl 


where 


rl  =  inner  radius  of  cell 
r2  =  outer  radius  of  cell 
zl  =  shorter  axial  distance 
z2  =  longer  axial  distance 

Analytically  integrating  with  respect  to  r  we  obtain  the  expression 

|/0)(r(  r2_)_r(yl_ 

fi  2  V2g(z)  V2 g(z) 

^(r22  -rl2)(z2-  zl) 


))dz 


(3) 


where 


T(x)  =  the  gamma  function 

Numerically  integrating  equation  (3)  using  a  step  integration  method  we  obtain  the 
expression 

(r(  r— )  -  T(  rl  )) 

_  2  V2g(zl)  V2g(zl)  ( 

—  (r22  -  rl2) 

2 
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for  the  plasma  density  in  a  cell  with  length  ( zl  -  z\ )  and  width  (r2-rl ).  Equation  (4)  is 
used  to  generate  the  plasma  density  tables  used  by  the  EWMPM  program,  which  in  turn 
allows  the  EMWPM  program  to  determine  the  refraction  and  attenuation  of  an  EM  wave 
as  a  result  of  the  plasma  distribution  calculated  by  the  EBS  program.  A  contour  plot  of 
the  plasma  density  resulting  from  equation  (4)  is  shown  in  Figure  34.  The  contour  plot 
represents  a  two  dimensional  slice  of  the  plasma  density  generated  by  the  electron  beam 
generator.  The  Gaussian  shape  of  the  electron  distribution  in  the  radial  direction  and  the 
natural  dispersion  due  to  electrons  at  larger  radial  distances  from  the  center  of  the 
electron  beam  being  distributed  over  a  larger  volume,  results  in  a  plasma  that  is  very 
dense  near  the  center  of  the  beam  and  decreases  in  density  very  rapidly  in  the  radial 

I 

direction.  The  shape  of  the  density  contour  plot  is  a  function  of  the  electron  energy  and 
not  the  electron  beam  current,  however,  the  number  of  electrons  and  hence  the  electron 
density  is  a  function  of  the  electron  beam  current.  Therefore,  the  plasma  density  scales 
linearly  with  electron  beam  power  for  electrons  at  a  constant  energy.  Also  note  that  the 
STD  of  the  electron  beam  achieves  a  maximum  value  at  an  intermediate  range. 

Similarly,  the  longitudinal  profile  achieves  a  maximum  in  the  vicinity  of  this  same  range. 
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a) 


b) 


Figure  34.  Contour  Plot  of  Electron  Beam  Generated  Plasma 
a)  500  keV  Electron  Beam  b)  300  keV  Electron  Beam 


EM  Wave  Attenuation  and  Refraction  due  to  the  Plasma  Density  Distribution 

This  section  presents  the  results  obtained  from  the  EMWPM  program  on  the 
refraction  and  attenuation  of  an  EM  wave  by  an  electron  beam  generated  plasma  with  and 
without  plasma  chemistry  losses.  The  refraction  of  an  EM  wave  by  an  electron  beam 
generated  plasma  is  presented  in  the  first  subsection.  The  last  two  subsections  contained 
in  this  section  summarize  the  average  attenuation  achieved  by  the  plasma  with  or  without 
losses  over  a  range  of  EM  wave  frequencies. 


Distortion  of  the  Electromagnetic  Wave  Front 

The  purpose  of  this  section  is  to  demonstrate  the  capability  of  the  EMWPM 
program  to  determine  the  refraction  of  an  EM  wave  with  finite  spatial  extent  as  it 
traverses  an  electron  beam  generated  plasma.  The  analysis  is  limited  to  a  plasma 
generated  by  a  single  electron  beam  at  1  MeV  electron  energy  and  an  air  density 
corresponding  to  an  altitude  of  5  km. 
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The  main  purpose  of  refraction  in  the  EMWPM  model  is  to  assure  the  path  of  the 
ray  through  the  plasma  is  accurate  allowing  for  accurate  calculations  of  the  plasma 
frequency.  The  higher  accuracy  of  the  plasma  frequency  along  the  path  EM  wave,  results 
in  a  higher  accuracy  in  the  attenuation  calculations  for  the  plasma.  The  analysis  provided 
below  is  a  qualitative  assessment  of  the  results  of  the  EMWPM  program. 

The  index  of  refraction  of  the  plasma  is  highly  dependent  on  the  EM  wave  and 
plasma  frequencies,  and  to  some  degree  the  collision  frequency  of  the  plasma  (Chapter  I. 
equation  (11)).  In  general,  the  closer  the  plasma  frequency  is  to  the  EM  wave  frequency, 
the  lower  the  index  of  refraction  of  the  plasma.  The  lower  the  index  of  refraction  and/or 

the  larger  Vn ,  the  more  that  the  wave  will  be  refracted  in  the  direction  of  the  gradient  of 
the  plasma  density  (see  Chapter  II,  equation  (1)),  which  in  general  is  pointing  away  from 
the  electron  beam  source  for  the  case  of  the  electron  beam  generated  plasma.  This 

happens  because  the  radius  of  curvature  of  the  EM  wave  is  proportional  to  Vn  and 
inversely  proportional  to  n.  If  the  EM  wave  frequency  is  less  than  the  plasma  frequency 
then  the  EM  wave  will  be  reflected.  If  the  EM  wave  is  reflected  then  the  plasma  will  not 
attenuate  the  EM  wave  because  it  does  not  propagate  through  the  plasma. 

In  Figure  36,  the  attenuation  of  an  EM  wave  over  the  radial  range  of  0  to  5  m  at 
various  frequencies  is  shown.  The  refraction  analysis  is  only  done  over  half  the  plasma, 
because  the  plasma  distribution  is  symmetric  about  the  x  axis.  In  Figure  37.a,  an  EM 
wave  at  a  frequency  of  300  MHz  traverses  a  plasma  predicted  by  the  EBS  program.  The 
plasma  frequency  varies  between  400  MHz  to  7  GHz  in  the  radial  range  between  0  and  5 
m,  therefore,  the  300  MHz  EM  wave  reflects  off  the  plasma.  For  the  1  GHz  frequency 
shown  in  Figure  36. b,  the  EM  wave  propagates  through  edges  of  the  plasma  unrefracted 
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where  the  plasma  frequency  is  lower.  As  the  plasma  frequency  increases  towards  the 
center  of  the  plasma,  the  refraction  of  the  EM  wave  is  greater.  At  the  center  of  the 
plasma,  the  EM  wave  is  reflected  because  the  plasma  frequency  is  greater  than  the  EM 
wave  frequency.  At  5  GHz,  the  EM  wave  penetrates  the  majority  of  the  plasma,  but  is 
noticeably  refracted.  As  the  frequency  of  the  wave  becomes  larger  than  5  GHz,  the  EM 
wave  is  refracted  even  less  resulting  in  very  little  distortion  of  the  EM  wave  front. 
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—  ray  Beam  Axis  -  x  (cm) 

Figure  35.  Refraction  in  the  Electron  Beam  Plasma  at  5  km  Altitude 
100  kW,  and  Electron  Energy  of  1  MeV  for  Frequencies  of  a) 

300  MHz  b)  1  GHz  c)  5  GHz  d)  10  GHz  e)  20  GHz  f)  30  GHz 

Spatial  Variations  in  Attenuation 

The  purpose  of  this  section  is  to  demonstrate  the  capability  of  the  EMWPM 
program  to  quantify  the  power  attenuation  that  occurs  over  the  width  of  the  incident 
wave.  The  analysis  is  limited  to  a  plasma  generated  by  a  single  electron  beam  with 
energies  between  300  keV  and  1  MeV  at  5  km  altitude.  A  more  in  depth  study  should  be 
performed  to  better  characterize  the  attenuation  properties  of  the  plasma  in  different 
environments. 

The  plasma  density  of  the  electron  beam  generated  plasma  decreases  as  e  in 
the  radial  direction,  which  results  in  a  large  variation  in  the  plasma  density  over  a  short 
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radial  distance.  The  amount  of  power  attenuation  over  the  width  of  the  incident  wave 
front  at  various  frequencies  is  shown  in  Figure  37.  The  graphs  in  Figures  35  and  37  give 
the  complete  picture  of  the  refraction  and  attenuation  resulting  from  a  1  MeV  electron 
beam.  In  Figure  37.a,  the  attenuation  for  a  wave  frequency  of  300  MHz  is  shown.  The 
amount  of  attenuation  is  quite  low  because  the  EM  wave  was  reflected  due  to  the 
frequency  being  lower  than  the  plasma  frequency.  In  Figure  37.b,  low  attenuation  is 
observed  within  a  100  cm  of  the  center  of  the  electron  beam  because  the  frequency  was 
less  than  the  plasma  frequency  in  that  range  as  well.  However  at  125  cm  and  greater 
radial  distance,  the  rays  are  refracted  rather  than  reflected,  therefore  they  experience 
between  -7  to  -24  dB  attenuation.  The  rays  past  125  cm  are  refracted  such  that  they 
propagate  through  a  lower  density  region  of  the  plasma  and  hence  attenuate  less.  In 
Figure  37. c  no  rays  are  reflected  immediately,  therefore  all  rays  attenuate.  The  center 
rays  are  attenuated  the  most  because  the  plasma  density  is  the  highest  near  the  center  and 
falls  of  exponentially  towards  the  edges  of  the  plasma.  The  EM  wave  in  Figure  37  .c  is 
attenuated  more  than  any  other  frequency  because  it  is  closest  to  the  plasma  frequency. 
This  relationship  is  clearly  seen  in  Figure  36,  where  the  complex  index  of  refraction 
reaches  a  maximum  value  when  the  ratio  of  the  EM  wave  frequency  to  the  plasma 

frequency,  0)1 0)p  ,  is  approximately  equal  to  one.  Also  Figure  36  indicates  that  for  the 

ratio  of  CO  I  CD p  greater  than  1.5  the  complex  index  of  refraction  asymptotically 

approaches  zero,  hence  resulting  in  negligible  attenuation  at  EM  wave  frequencies  much 
greater  than  the  plasma  frequency. 
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Figure  36.  Plot  of  the  Imaginary  Index  of  Refraction, 
n„  versus  0)1  (0p  with  a  vl (Op  of  0.5 

For  frequencies  in  excess  of  5  GHz,  the  attenuation  decreases  steadily  as  the  EM 
wave  frequency  becomes  much  larger  than  the  plasma  frequency,  which  corresponds  with 
the  behavior  of  the  imaginary  index  of  refraction  shown  in  Figure  36. 
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Figure  37.  Spatial  Attenuation  at  5  km  Altitude 

for  Wave  Frequencies  of  a)  300  MHz  b)  lGFb  c)  5  GHz 

d)  10  GHz  e)  20  GHz  f)  30  GHz 


The  analysis  of  the  spatial  attenuation  of  the  EM  wave  was  performed  for  the 
initial  electron  energies  between  300  keV  and  1  MeV,  at  a  fixed  electron  beam  power,  to 
gain  some  insight  into  the  attenuation  resulting  from  different  density  profiles  of  the 
plasma.  For  the  EM  wave  frequency  analysis,  the  spatial  attenuation  of  the  plasma  was 
averaged  over  a  radius  of  5  m  to  obtain  an  average  attenuation  factor.  The  attenuation 
factor  was  then  plotted  versus  frequency  for  several  electron  beam  energies.  Figure  38 
indicates  that  attenuation  is  larger  at  lower  EM  wave  frequencies  when  the  electron  beam 
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has  higher  electron  energies  because  the  plasma  density  is  lower.  For  higher  frequencies, 
lower  electron  energies  are  needed  for  greater  attenuation,  because  plasma  density 
increases  due  to  the  smaller  plasma  volume,  hence  the  plasma  frequencies  are  nearer  to 
the  higher  EM  wave  frequencies.  However,  there  is  a  point  where  lower  electron 
energies  result  in  a  beam  with  a  restrictive  radial  extent.  If  attenuation  at  higher 
frequencies  is  desired  then  the  electron  beam  generator  may  be  operated  at  lower  electron 
energies,  resulting  in  only  a  portion  of  the  EM  wave  near  the  centerline  being  attenuated 
and  as  a  result  any  attenuation  in  the  radial  wings  is  negligible.  This  solution  results  in 
larger  attenuation  of  the  EM  wave  over  a  broader  band  of  frequencies.  Whereas  the 
higher  electron  energies  provide  a  broader  spatial  coverage,  but  the  range  of  highly 
attenuated  frequencies  is  limited. 
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Figure  38.  Attenuation  versus  Frequency  at  5  km  Altitude 
for  Initial  Electron  Energies  of  a)  1  MeV  b)  750  keV 
c)  500  keV  d)  300  keY 
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Spatial  Variations  in  Attenuation  Considering  Plasma  Loss  Mechanisms 

In  the  Plasma  Density  Loss  Mechanisms  section  of  Chapter  IV,  it  was  shown  that  the 
plasma  reached  an  approximate  steady  state  in  a  few  microseconds  at  all  electron  and  ion 
density  injection  rates  of  interest.  As  a  result,  the  effects  of  electron  attachment,  detachment, 
and  recombination  on  the  plasma  density  were  approximated  as  a  constant  loss  factor  to  the 
plasma  density  referred  to  as  the  Plasma  Density  Loss  Factor  (PDLF).  Therefore,  to  determine 
the  plasma  density  distribution  with  loss  mechanisms,  we  can  simply  multiply  the  plasma 
density  without  loss  mechanisms  by  the  PDLF.  From  Figure  29,  we  see  that  the  lowest  PDLF 
is  achieved  by  operating  the  electron  beam  for  0.1  ms,  which  results  in  a  PDLF  of  0.00134. 
Since  in  the  previous  section  the  electron  beam  was  operated  for  0.0005  s  at  100  kW,  we  must 
increase  the  power  by  a  factor  of  5  to  maintain  the  same  pulse  energy  resulting  in  an  electron 
beam  power  of  500  kW.  The  plasma  density  distribution  resulting  from  the  application  of  the 
attenuation  factor  results  in  a  negligible  amount  of  attenuation  (less  than  10  dB  for  all 
frequencies). 

Summary  and  Conclusions 

The  objective  of  this  study  was  to  develop  a  suite  of  computational  tools  to 
analyze  the  propagation  of  an  electromagnetic  wave  in  a  plasma  generated  by  injecting  a 
relativistic  electron  beam  into  the  air.  The  suite  had  three  major  physical  components: 
electron  beam  propagation  and  plasma  generation,  evolution  of  the  plasma  densities  due 
to  plasma  chemistry,  and  EM  wave  propagation.  These  components  were  translated  into 
three  separate  computer  programs.  Data  generated  in  the  programs  were  post-processed 
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using  Mathematica,  which  provided  analytic  fits  to  source  terms,  particle  densities,  and  a 
graphical  interpretation  of  system  characteristics. 

The  spatial  extent  and  density  distribution  of  the  thermal  plasma  source  generated 
by  injecting  a  relativistic  electron  beam  into  the  air  were  investigated  using  a  stochastic 
approach  based  on  an  axisymmetric  Monte  Carlo  model.  Establishing  a  reliable  set  of 
cross  sections  for  this  simulation  was  crucial.  A  thorough  review  of  the  literature 
resulted  in  an  experimentally  consistent  set  appropriate  for  this  analysis:  species  specific 
and  relativistically  correct.  The  differential  and  total  cross  sections  for  scattering  and 
ionization  of  Mott,  Bethe,  and  Kim  were  analyzed  and  incorporated  into  the  Monte  Carlo 
simulation.  The  Monte  Carlo  calculation  was  validated  using  a  limited  case  in  which 
electrons  only  experienced  forward  scattering  and  were  limited  to  losing  an  average 
ionization  energy  per  ionization  process.  The  results  of  a  limited  study  were  consistent 
with  independent  analytic  and  numerical  implementations.  Once  the  Monte  Carlo 
simulation  was  validated,  the  data  associated  with  the  general  solution  was  analyzed  and 
the  parametric  dependence  of  the  source  footprint  explored.  A  limited  exploration  of  the 
dependence  of  the  plasma  distribution  on  neutral  densities  and  the  initial  electron 
energies  provided  scaling  relations.  The  longitudinal  extent  of  the  plasma  varied  linearly 
with  respect  to  the  energy  of  the  electron  beam,  as  did  the  transverse  extent.  Both  the 
longitudinal  and  transverse  extent  varied  inversely  with  respect  to  the  neutral  density, 
which  translates  into  an  altitude  dependence  if  an  exponential  atmosphere  model  is 
employed.  For  a  density  associated  with  an  altitude  of  5  km,  the  plasma  longitudinal 
extent  ranged  from  52  to  868  cm  and  the  standard  deviation  of  the  transverse  radial  extent 
ranged  from  18  to  292  cm  for  initial  electron  energies  between  100  keV  and  1  MeV. 
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The  plasma  source  distribution  from  the  Monte  Carlo  model  was  used  in  the 
second  code,  a  spatial  and  temporal  solution  of  the  plasma  evolution  based  on  a  time 
dependent  analysis  of  the  plasma  rate  equations  in  air.  Only  volumetric  processes  were 
considered.  A  set  of  relevant  kinetic  rate  coefficients  and  species  rate  equations  were 
assembled  and  a  time-dependent,  spatially  resolved  solution  of  the  plasma  densities  was 
achieved  using  a  numerical  integration  of  the  coupled  equations.  For  the  conditions 
examined,  a  pseudo-steady  state  of  the  electron  density  is  achieved.  The  resulting  plasma 
is  a  three  component  plasma  consisting  of  electrons,  positive  ions  and  negative  ions.  The 
negative  ions  of  molecular  oxygen  are  the  majority  negative  species  due  to  the  rapid 
three-body  attachment  process  at  the  neutral  densities  considered.  The  rate  equation 
results  were  interpreted  by  employing  a  simplified,  steady  state  analysis.  The  values  of 
the  electron  and  negative  ion  density  were  established  analytically  and  the  scaling  with 
respect  to  electron  beam  current  explored  and  explained.  For  the  range  of  electron  beam 
pulse  lengths  between  5  ms  and  0.1  ms,  the  density  of  the  electrons  in  an  air  plasma 
produced  by  an  electron  beam  will  be  reduced  by  a  factor  between  3x10'  to  2x10' 
respectively,  due  to  attachment  and  recombination  processes. 

The  third  code  developed  evaluates  the  attenuation  and  refraction  of  an  EM  wave 
in  a  plasma  of  arbitrary  spatial  distribution.  A  ray  tracing  method  based  on  the  eikonal 
approach  of  Sommerfeld  was  implemented  numerically  and  validated  against  analytic 
solutions  relating  to  radio  wave  propagation  in  the  ionosphere.  This  approach  enabled 
evaluation  of  both  wave  refraction  and  absorption.  Neglecting  attachment,  the  resulting 
plasma  was  also  found  to  significantly  refract  and  attenuate  the  EM  waves  at  reasonable 
electron  beam  power  settings  and  pulse  lengths.  The  amount  of  refraction  was  very 


137 


dependent  on  the  proximity  of  the  plasma  and  the  EM  wave  frequencies.  If  the  plasma 
frequency  over  the  transverse  extent  of  the  plasma  was  much  larger  or  smaller  than  the 
EM  wave  frequency  then  very  little  refraction  occurred.  However,  if  the  values  of  the 
plasma  and  EM  wave  frequencies  were  close,  then  the  EM  wave  was  refracted  by  a 
significant  amount. 

In  summary,  this  study  successfully  integrated  plasma  generation,  plasma 
evolution  and  wave  propagation  analyses  to  permit  a  quantitative  evaluation  of  the  effects 
of  an  electron  beam  generated  plasma  on  an  EM  wave.  The  EBS  simulation  was  capable 
of  characterizing  the  plasma  density  distribution  resulting  from  a  relativistic  electron 
beam.  The  EMWPM  program  was  also  capable  of  determining  refraction  and  attenuation 
of  an  electromagnetic  wave  traversing  an  arbitrary,  collisional  plasma.  The  plasma 
chemistry  model  also  demonstrated  its  ability  to  analyze  the  temporal  evolution  of  the 
plasma  due  to  the  chemical  processes  in  the  plasma  and  translate  those  results  into  a 
reduction  of  the  plasma  density  distribution  over  time.  These  models  are  very  flexible 
and  should  be  used  to  further  examine  the  parameter  space. 

Limitations  of  the  Study  and  Recommendations 

For  the  plasma  chemistry  calculation  the  change  in  energy  of  the  thermal 
electrons  from  high  to  low  energies  due  to  inelastic  collisions  with  molecules  and  ions 
was  modeled  very  coarsely  and  a  more  refined  calculation  should  be  used  such  as  a 
Boltzmann  transport  calculation.  The  electron  temperature  calculations  are  important 
because  the  rate  constants  for  attachment  and  detachment  processes  are  very  dependent 
on  the  temperature  of  the  electrons  (such  as  attachment  reactions  3  and  4  in  Appendix  C). 
Given  the  importance  of  attachment  loss  mechanisms  and  the  dependence  of  these 
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processes  on  the  energy  distribution  of  the  thermal  electrons,  a  refined  treatment  of 
thermal  electrons  is  appropriate.  Spatial  variations  in  the  electron  temperature  through 
out  the  plasma  were  also  not  considered  in  the  study.  Spatial  variations  in  the  electron 
temperature  will  result  in  different  collision  frequencies  and  hence  a  different  amount  of 
attenuation  of  the  EM  wave.  Any  heating  of  the  air  by  the  electron  beam  was  also 
considered  negligible.  However,  this  may  not  be  true  if  the  electron  beam  is  operated  at  a 
higher  power  and  for  a  longer  duration  than  was  considered  in  this  study.  Also,  the  only 
constituent’s  of  the  air  considered  in  the  rate  calculations  of  the  plasma  density  were 
oxygen  and  nitrogen.  Other  molecules,  even  though  they  are  found  in  small 
concentrations,  may  have  a  significant  impact  on  the  densities  of  the  plasma  constituents. 

Experiments  to  determine  the  density,  spatial  distribution,  and  temporal  evolution 
of  the  plasma  should  be  performed  to  verify  the  results  of  the  EBS  and  plasma  chemistry 
programs.  Measurements  of  the  power  attenuation  due  to  a  collisional  plasma  could  also 
be  performed  to  verify  the  EMWPM  program  and  assess  wave  propagation  in  a  general 
class  of  artificially  generated  plasmas. 

Other  recommendations  for  future  work  include:  A  thorough  investigation  of  the 
parameter  space  of  the  electron  beam  densities  and  the  resulting  attenuation  and 
refraction.  Also  the  effects  on  the  complex  index  of  refraction  due  to  a  three-component 
plasma  should  be  investigated  due  to  the  majority  of  the  negatively  charged  species  in  the 
air  plasma  being  molecular  oxygen  ions. 
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Appendix  A:  Derivation  of  the  Radius  of  Curvature 


Starting  with  the  scalar  wave  equation 


V2w  +  k2u  —  0 


and  assuming  that  the  solution  to  equation  (1)  is  given  by 


u  =  Ae  ik»s 


K  —  ^£0Mo  ® 


where 


A  =  amplitude  factor 
S  =  the  eikonal  function 

where  we  consider  u  to  be  a  rapidly  varying  function  of  position  and  A  and  S  as  slowly 

varying  functions  of  position  (Sommerfeld,  1964:330).  Substituting  equation  (2)  into  the 

wave  equation  we  obtain 
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where  the  terms  indicated  by . do  not  become  infinite  as  ko  approaches  infinity. 

Equation  (1)  is  satisfied  if  S  and  A  satisfy  the  differential  equations 


'XT' 


+ 


fasY  fasY 
Kdyj 


- n 2  =0,  n  = 


j V2S+VLogA  VS  =0 


(4) 

(5) 


n  =  index  of  refraction 

Equation  (4)  is  referred  to  as  the  differential  equation  of  the  eikonal.  Equation  (5)  does 
not  require  that  the  gradient  of  the  logarithm  of  A  be  perpendicular  to  the  gradient  of  S; 
therefore,  it  permits  discontinuities  of  A  in  these  directions.  If  S  is  equal  to  a  constant 

value,  such  as  x2  +  y2  +  z2  =  C2 ,  then  the  function  u  is  at  a  constant  phase,  which 
results  in  a  wave  surface.  The  normal  to  the  wave  surface,  the  gradient  of  S,  is  the  ray 
propagation  direction.  If  the  index  of  refraction  varies  in  space  then  the  rays  will  curve  in 
accordance  with  equation  (4). 

In  an  optically  homogenous  medium  the  simplest  solution  to  equation  (4)  is  the 
linear  function 

S  =  n(ax  +  /3y  +  yz)  (6) 

where 

1  =  a1  +  /32  +  Y2  (7) 

This  solution  indicates  the  waves  surfaces  are  planes,  the  propagation  direction  is 
perpendicular  to  the  plane  described  by  equation  (6),  and  the  rays  are  parallel  which  is  in 
accordance  with  geometric  optics.  Other  solutions  for  an  optically,  homogenous  medium 
include  spherical  and  cylindrical  wave  fronts. 
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The  unit  vector  in  the  direction  of  the  ray  propagation  is  given  by  the  expression 


-VS 

s  =  — — r  =  —VS 


VS 


n 


(8) 


The  curl  of  the  gradient  of  any  function  is  defined  as  zero.  This  fact  gives  us  the 
expression 


V  x  (VS)  =  V  x  (ns)  =  0  (9) 

According  to  Sommerfeld,  “This  condition  is  equivalent  to  the  existence  of  the  eikonal. 
All  ray  bundles  (straight  or  curvilinear)  realized  in  geometrical  optics  are  normals  to 
surfaces  and  are  distinguished  from  more  general  systems  of  curves  in  that  they  satisfy 
the  condition  (9).” 

From  Stokes’  theorem  we  obtain  the  integral  form  of  equation  (9) 


which  results  in 


V  x  (ns)  •  da  =  ^(n?)  •  ds  =  0  (10) 


p  2 

|  ns  ■  ds  =  S2  -  Sl  (11) 


pi 

This  indicates  that  the  change  in  the  eikonal  is  independent  of  the  path  of  the  ray. 

From  Figure  3,  we  have  the  incoming  ray,  s  ,  tangent  to  the  circle  at  point  P  and  we  have 
the  refracted  ray,  s' ,  at  P’.  The  curvature  of  the  ray  is  defined  as  the  angle  between 
s  and  s'  divided  by  the  distance  PP’ 


K 


1  _  0 

R  ~  PP' 


(12) 


Since  |?|  =  1 ,  the  angle,  0  ,  is  equal  to  s  -  s'  and  for  a  very  small  change  in  direction, 
we  can  make  the  statement  that  PP'=  ds  .  Therefore,  we  define  curvature  as 
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(13) 


ds  =  direction  of  the  curvature  vector 

Using  the  chain  rule  on  equation  (13)  in  Cartesian  coordinates,  we  obtain 

ds  _ds  dx  +  ds  dy  +  ds  dz 

ds  dx  ds  dy  ds  dz  ds 

dx  dy  dz  „  ,  _  _ 

where  the  full  derivative  — ,  — ,  —  are  the  components  of  the  s  vector.  Therefore, 

ds  ds  ds 

the  chain  rule  can  be  written  in  the  following  form 

ds  ds  ds  ds 

*"*'*  +  *'’  +  fc'*  (15) 

.  i  .  .2 

Since,  ?  =  1 ,  then  U'l  =1  and  by  taking  the  gradient  of  both  sides  of  the  equation 


^V|?|2  =  3xVsx  +  syWsy  +  szWsz  =  0 


subtracting  equations  (15)  and  (16)  results  in 

ds  .ds  -  .  ,  .ds  .  ,  .ds  - 

The  x-component  of  this  vector  equation  is  given  by 


ds  f  dsx  dsy  |  ( dsx  dsz 


which  using  vector  identities,  it  can  be  shown  that 


(Vxi)xi 


By  using  the  identity 


Vx(/A)  =  yVxA- AxV/ 
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with  the  fundamental  equation  (9)  which  characterizes  the  ray  vector,  s  ,  we  obtain  the 
expression 

Vx?  =  — sxVn  (21) 

n 

If  we  substitute  equation  (21)  into  (19),  we  obtain  the  equation 

K=—  (sxVn)xs  (22) 

n 

Using  the  “B AC  CAB”  triple  product  rule  we  obtain  the  final  expression  for  the  curvature 
vector  of  a  ray  in  an  inhomogeneous  medium 

K  =  —  (Vn  -  s(s  •  Vn))  (23) 

n 

From  this  we  see  that  the  principle  normal  K ,  the  tangent  s  ,  and  the  Vn  all  lie  in  one 
plane  (Sommerfeld,  1964:339). 
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Appendix  B:  Input  Parameters  for  the  EBS  Simulation 

This  appendix  provides  a  description  of  the  input  variables  and  files  used  in  the 

EBS  simulation.  The  input  file  format  for  the  EBS  program  is  the  standard  Fortran  90 

NAMELIST  I/O  format.  The  following  is  an  example  of  the  format  of  the  EBeam.inp  file 

!  Electron  Beam  Simulation  Input 
&start  Runs  =  1  , 

Power(l)  =  30000.0,  30000.0  , 

EEnergy(l)  =  100000.0,  1000000.0, 

Altitude(l)  =  10000.0,  10000.0, 

MaxCells(l)  =  200.0,  200.0,  200.0,  200.0,  200.0  , 

CellSize(l)  =  7.5, 7.5, 

SwapSize(l)  =  20000,  20000, 

MoveNewE(l)  =  .TRUE. ,  .TRUE. , . 

OutputFile(l)  =  ‘EBeamOutputl.dat’  ,  ‘EBeamOutput2.dat’  , 

NumSimE(l)  =  100,  100  / 

where 

!  -  Comment 

&  -  start  of  the  NAMELIST  input 
start  -  the  name  of  the  NAMELIST 

Runs  -  variable,  can  be  any  type  integer,  real,  logical,  character,  etc. 

Power(l)  -  Array,  can  be  an  integer,  real,  logical,  or  character  array.  (1)  indicates 
that  data  input  starts  in  the  first  element  of  the  array. 

/  -  end  of  the  NAMELIST  input 

The  input  files,  Default_File.dat  and  Molecule_Data_File.dat  for  the  EBS  program  are 
the  same  format,  however,  only  the  variables  that  can  be  used  in  the  particular  file  are 
different.  A  list  of  all  the  variables  used  in  the  EBeam.inp  file  are  shown  in  Table  14. 

The  first  column  contains  the  name  of  the  variable,  the  second  column  the  units  for  the 
variable,  a  brief  description  of  the  variables  is  included  in  the  third  column,  and  the 
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fourth  and  fifth  columns  contain  the  minimum  and  maximum  values  that  will  be  accepted 
by  the  EBS  simulation. 


Table  14.  Input  Variables  to  EBS  program 


Ebeam.inp  Run  Parameters 

Name 

Units 

Description 

Min 

Max 

Power 

W 

Electron  beam  (Ebeam)  power 

0.0 

10A9 

EEnergy 

eV 

Initial  energy  of  the  electrons 

0.0 

10A6 

AirNumDens 

#/cmA3 

Number  density  of  the  air 

0.0 

10A33 

Altitude 

m 

Altitude  of  the  EBeam 

0.0 

10A5 

BeamRadius* 

cm 

Radius  of  electron  beam  nozzle 

10A-3 

10A4 

PRI 

s 

Pulse  Repetition  Interval  (PRI)  of  the 
electron  beam 

0.0 

10 

DutyCycle 

# 

Percent  time  that  Ebeam  is  on  during  the 
pulse  interval 

0.0 

1.0 

Efficiency 

# 

The  efficiency  at  which  the  electron  gun 
converts  power  into  electrons 

0.0 

1.0 

NumSimE** 

# 

Number  of  initial  electrons  in  the 
simulation 

0 

10A5 

IonEnergy 

eV 

Average  ionization  energy  of  the  atoms  or 
molecules  in  the  simulation 

0.0 

10A4 

MaxCells** 

# 

Creates  a  grid  that  is  MaxCells  long  in  the 
axial  direction  (X)  by  2*MaxCells  in  the 
Transverse  Direction  (Y) 

0 

1000 

CellSize 

cm 

Size  of  one  side  of  a  cell  in  the  grid 

0.0 

10A3 

RunTime 

s 

Maximum  amount  of  time  that  the 
simulation  will  run  for  each  electron 

0.0 

0.01 

NumSteps* 

# 

Number  of  intermediate  parameter  values 
between  the  min  and  max  value  when  in 
analysis  mode 

0.0 

100.0 

NRGSTD* 

eV 

Standard  deviation  in  the  distribution  of 
the  electron’s  initial  energy  distribution 

0.0 

10A6 

AngleSTD* 

Radians 

Standard  deviation  in  the  distribution  of 
the  electron’s  initial  angle 

0.0 

Pi/2 

Logical 

If  true  all  scattering  events  are  Isotropic 

Logical 

If  true  simulates  Ionizing  events 

Logical 

If  true  simulates  Elastic  Scattering 

InitDist* 

Logical 

If  false  electron  beam  is  Monoenergetic,  and  all  electrons 
have  an  initial  angle  of  0.0  rad 

If  true  electron  beam  has  a  gaussian  distribution  using 
NRGSTD  and  AngleSTD 

146 


NoEjectNRG 

Logical 

If  true  ejected  electrons  have  no  energy  after  being  ionized 
from  the  molecule 

ExpAtm 

Logical 

If  true  use  the  exponential  atmosphere  model  to  determine 
air  number  density 

If  false  use  AirNumDens  value  for  air  number  density 
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Appendix  C:  Primary  Chemical  Reactions  in  the  Plasma 

Table  14  provides  a  list  of  all  the  chemical  kinetic  equations  that  were  used  in  the  plasma 
chemistry  calculations.  Column  1  is  the  reaction  number.  Column  2  is  the  reaction 
process  itself.  Columns  3  and  4  provide  the  minimum  and  maximum  values  of  the  rate 
constant  for  the  molecular  temperature  or  average  electron  energy  range  listed  in  column 
5.  Column  6  is  the  source  of  the  data  on  the  rate  constant. 

Table  15.  List  of  Dominant  Chemical  Reactions  for  a  Nitrogen-Oxygen  Plasma 


# 

Reaction  Process 

Rate 

Constant  or 
Cross- 
Section 

Rate  Constant 
or  Cross- 
Section 

Temp, 
or  Avg 
e' 

Energy 

Source 

Attachment  -  Two  Body 

Min 

Max 

1 

c  +  02  — ^  O  0 

0 

1.25xl01K 

cm 

at  6.7  eV 

3.6  to 

12  eV 

McDaniel 

p410 

1 

6  +  02  — >  02 

2xl0'19 

Niles, 
Table  V, 
React.  17 

2 

e  +  0^>0~ 

■ 

1 

Niles 
React.  15 

■ 

e  +  N2-^N  +  N~ 

Nil 

Nil 

Alpert 
p  99 

■ 

e  +  N  — >  N  ~ 

Nil 

Nil 

■a 

Attachment  -  Three  Body 

3 

e  +  20 2  — ^  02  +  02 

HI 

McDaniel 

p411 

■ 

195-600 

K 

Niles 

e  +  02  +  N  2  — >  N  2  +  02 

lxlO132 

0.1  to  1 
eV 

McDaniel 
p  409 

lxKT*1 

300K 

Niles 
Table  V, 

5 

e  +  02  +  0  — )  0  +  02 

IxitF 

Niles 

J.  Chem 
Phys 
52:408 
(1970) 
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Electron  Recombination 
Two  Body 


e  +  N2  — »  N2 


46  e  +  N^_^2N 

Electron  Recombination  - 
_ Three  Body _ 

10  e  +  02+  +M  ^02+M 


->0  + 


Electron  Detachment  - 
Two  Bod- 


O  +0  — ^  C)2  "I-  € 


02  +  O  — >  Oj  +  e 


O  +  02  — ^  Oj  +  c 


02  +  02  *  — ^  202  +  e 


1x10 


2.1x10' 


1.2x10' 


Dettmer 
p  207 


Niles, 
Table  V, 
React.  11 
Niles, 
Dettmer 
React.  19 
Alpert 
99 


pert 

99 


Niles 


Niles 

J.  Chem 
Phys 
52:408 
(1970) 

Niles 

J.  Chem 
Phys 
52:408 
(1970) 

Dettmer 
p  263 
React.  3 


Dettmer 
p  264 

16 

Niles 
React.  9 

Dettmer 
p  270 
React.  117 

Dettmer 
p  270 
React.  118 
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02  "t"  02  — ^  202  +  £ 


O  +  2VO  — »  1V02  +  e 


Additional  02  Losses 


02  +  O  — ^  O  +  02 


17  02++0_+M^03+M  2x10 


02  ■+■  02  — ^  02  +  20 


O2  +  02  — ^  20j 


02  +  N2  ->  NO+  +  NO 


02  +N^NO+  +0 


Additional  02+  Gains 
0+  +  0  +  M  ^02+  +M 


21  o++o2->o2++o 


02  +  N2+  — >  02+  +  N2 


2.0x10 


Niles 
React.  11 

69:81 

Dettmer  p. 
263 

React.  9 

69:82 

Dettmer  p. 
263 

React.  12 

70 

Dettmer  p. 
264 

React.  29 

69:81 

Dettmer  p. 
265 

React.  43 

Niles 

React.  119 

Niles 

70  Dettmer 
p  264 
React.  30 


69:44  Dettmer 
p  265 

_ React.  39 

See 

Additional 

n2+ 

Losses 
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Additional  02  Losses 

(?2  "I"  O2  — ^  ^2  20 

See 

Additional 

<v 

Losses 

22 

0+  +  O2  — )  02  +  0 

l.OxlO'7 

69:81 

Dettmer 
p  264 
React.  23 

23 

0  +  02  — ^  02  +  0 

3.3xlO,u 

69:69 

Dettmer 
p  264 
React.  26 

02  +  02  — )  202 

— 

See 

Additional 

<v 

Losses 

Additional  02  Gains 

24 

0  +  02  — ^  02  ^ 

78:124 

7 

Dettmer 
p  270 
React.  114 

Additional  0+  Losses 

25 

0+  +0~  ->20 

2.7xl0‘v 

69:82 

Dettmer 
p  263 
React.  8 

26 

0+  ~\~  0  +0  — ^  02  ^  ^ 

2.0x10-*’ 

69:82 

Dettmer 
p  263 
React.  11 

0 +  +  02  — >  O  +  02 

See 

Additional 

<v 

Losses 

0+  +O  +  M  -^02  +0 

See 

Additional 
02  Gains 

Additional  0~  Losses 

0+  +0~  -*20 

See 

Additional 

0+ 

Losses 

02  +  ->  0  +  o2 

See 

Additional 

Losses 
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1 

0 +  +  0  +0  — >  02  +  0 

See 

Additional 

o+ 

Losses 

1 

$2  +  0  +  M  — ^  03  +  M 

See 

Additional 

Losses 

1 

0  +  02  +  M  — >  02  +  0 

See 

Additional 
02  Gains 

33 

o~  +n2+  ->n2+o 

2.0xl0'; 

Niles 

35 

0~  +  N2+  +  N2  ->  2N2  +  0 

Hi 

Niles 

Additional  0~  Gains 

1 

0  +  02  — ^  02  +  0 

See 

Additional 

<V 

Losses 

Additional  O  Losses 

0  +  e—>0++  2e 

See 

Additional 
0+  Gains 

20  +  02  — >  20 2 

69:93 

Dettmer 
p  264 
React.  28 

1 

0+  +  O  +  O  — )  02  +  0 

See 

Additional 

0+ 

Losses 

28 

20  +  M  —>  02  +  M 

7.4xl0’a3 

69:93 

Dettmer 
p  265 
React.  3 1 

29 

0  +  202  — )  O3  +  02 

6.0xl0"4 

69:93 

Dettmer 
p  265 
React.  32 

30 

20  +  02  — )  O3  +  0 

6.0xl0"4 

69:93 

Dettmer 
p  265 
React.  33 

152 


Additional  O  Gains 


02  +  O  — >  O  +  02 


02  +  02  — )  02  +  20 


0++02+M->02  +0 


02  +  02  — ^  02  +  20 


0+  +  02  — >  02  +  o 


O  +  02  — >  02  +  o 


0+  +0'  ^20 


0+  +02'  ->0  +  0 


See 

Additional 

<v 

Losses 

See 

Additional 

Losses 

See 

Additional 
02  Gains 

See 

Additional 

02+ 

Losses 

See 

Additional 

0; 

Losses 

See 

Additional 
02  Gains 

See 

Additional 

0“ 

Losses 

See 

Additional 

<v 

Losses 

See 

Additional 
0+  Gains 


See 

Additional 


38 

02  +  N2  ~^NO+  +NO 

Niles 

React.  123 

mm 

Niles 

React.  176 

40 

n  +  n2+  -*n2+n+ 

lxlO'12 

Niles 
React.  90 

41 

N2+  +0^  NO+  +  N 

42 

02  +  N 2  +  N2  — ^  02  +  2N 2 

2.0xl0'25 

Niles 
React.  66 
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Appendix  D:  Results  of  Rate  Equation  Calculations 


Figure  39  shows  a  sample  of  the  results  of  the  chemical  kinetics  calculations  with 
different  source  term  values.  The  graphs  in  Figure  39  represent  the  densities  of  the 
constituents  of  the  plasma  as  a  function  of  time.  The  graph  labeled  Ne  is  the  electron 
densities  in  the  low  thermal  electron  energy  range  (<  0.1  eV)  as  defined  in  Table  11.  The 
graphs  labeled  e*  and  e**  are  the  electron  densities  in  the  high  and  medium  energy 
ranges  respectively.  All  references  to  reaction  numbers  in  this  appendix  are  referring  to 
the  reactions  listed  in  Table  15.  For  the  calculations  presented  in  this  appendix  the 
electron  beam  was  turned  on  for  0.005  s  and  then  turned  off.  Therefore,  these  graphs 
represent  the  response  of  the  air  to  an  influx  of  electrons  and  ions  and  then  the  relaxation 
of  the  constituents  of  the  plasma  after  the  source  has  been  removed.  Figure  39  shows 
that  the  most  abundant  constituents  other  than  molecular  nitrogen  and  oxygen  are  the 
positive  and  negative  ions  of  molecular  oxygen.  The  negative  ions  of  molecular  oxygen 
have  a  high  concentration  primarily  because  of  the  three-body  attachment  processes  of 
electrons  with  molecular  oxygen  (reactions  3,  and  4)  in  the  low  and  medium  energy 
ranges.  The  positive  ions  of  molecular  oxygen  have  a  high  concentration  because  of 
charge  exchange,  reaction  37,  which  results  in  a  rapid  transfer  of  electrons  from  neutral 
oxygen  molecules  to  positive  nitrogen  ions  forming  positive  oxygen  ions  and  neutral 
nitrogen  molecules.  Atomic  oxygen  and  its  negative  ion  are  also  in  plentiful  supply 
because  they  are  a  product  of  the  very  fast  dissociative  attachment  process  given  by 
reaction  1.  This  analysis  is  true  for  the  range  of  source  terms  shown  in  Figure  39.  Other 
reaction  processes  only  start  becoming  dominant  at  higher  electron  beam  currents. 
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Figure  39.  Chemical  Kinetics  Results  for  Electron  Beam  Pulse  of  5  ms 
with  Source  Terms,  y ,  of  a)  2xl013  b)  2xl014  c)  2xl015 
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electron  beam  generated  plasma  in  refracting  and  attenuating  an  EM  wave.  The  spatial  extent  and  density  distribution  of  a  plasma  generated  by  a  relativistic  electron 
beam  was  determined  using  an  axisymmetric  Monte  Carlo  model.  This  plasma  density  distribution  was  used  as  a  source  term  in  the  second  code,  a  temporal  solution 
of  the  plasma  evolution  based  on  a  time  dependent  analysis  of  the  plasma  rate  equations.  The  third  code  developed,  evaluates  the  attenuation  and  refraction  of  an  EM 
wave  in  the  resulting  plasma  by  using  a  ray  tracing  method  based  on  the  eikonal  approach  of  Sommerfeld.  The  theoretical  foundation  and  validation  procedures  are 
presented  for  each  program.  A  limited  exploration  of  the  dependence  of  the  plasma  distribution  on  neutral  density  and  the  electron  beam  energies  was  performed.  At  a 
neutral  density  corresponding  to  5  km  of  altitude,  the  plasma  longitudinal  extent  ranged  from  52  to  868  cm  and  the  radial  extent  ranged  from  18  to  292  cm  for  initial 
electron  energies  between  100  keV  and  1  MeV  respectively.  Plasma  chemistry  plays  a  critical  role  in  determining  the  electron  plasma  density  and  dictates  the  beam 
format  required  to  achieve  a  desired  level  of  EM  wave  attenuation. 
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